Finite element method for radially symmetric solution of a multidimensional semilinear heat equation

TL;DR

This paper develops and analyzes a finite element method for solving radially symmetric solutions of multidimensional semilinear heat equations, providing error estimates and validating results through numerical examples.

Contribution

It introduces optimal error estimates for finite element solutions of semilinear heat equations with radial symmetry, including both symmetric and nonsymmetric formulations.

Findings

Optimal order error estimates in $L^ Infty$ and weighted $L^2$ norms.

Validation of theoretical results through numerical examples.

Abstract

This study aims to present the error and numerical blow up analyses of a finite element method for computing the radially symmetric solutions of semilinear heat equations. In particular, this study establishes optimal order error estimates in $L^\infty$ and weighted $L^2$ norms for the symmetric and nonsymmetric formulation, respectively. Some numerical examples are presented to validate the obtained theoretical results.

Figures

Equations

$$u_t=u_{xx}+\frac{N-1}{x}{u}_x+f(u),$$

$$u_x(0,t)=u(1,t)=0,$$

$$u(x,0)=u^0(x),$$

$$N\ge 2\text{mbox}integer$$

$$U_t=\Delta U+f(U),$$

$$U=0,$$

$$U(0,\mathbf{x})=U^0(\mathbf{x}),$$

$$U_t=\Delta U+U^{1+\alpha },\alpha >0$$

$$\forall \mu >0,\exists M_{\mu }>0:|f(s)-f(s^{\prime })|\le M_{\mu }|s-s^{\prime }|(s,s^{\prime }\in \mathbb{R},|s|,|s^{\prime }|\le \mu ).$$

$${\int }_I{x}^{N-1}{u}_t\chi dx+{\int }_I{x}^{N-1}{u}_x{\chi }_{x}dx={\int }_I{x}^{N-1}f(u)\chi dx.$$

$${\int }_{I}x{u}_t\chi dx+{\int }_I[x{u}_x{\chi }_x+(2-N)u_x\chi ]dx={\int }_{I}xf(u)\chi dx.$$

$$0=x_0<x_1<\cdots <x_{j-1}<x_j<\cdots <x_{m-1}<x_m=1,$$

$$S_h=\{v\in H^1(I)\mid v\in {\mathcal{P}}_1(I_j)(j=1,\cdots ,m),v(1)=0\}.$$

$${\varphi }_j(x_i)={\delta }_{ij},$$

$$t_0=0,t_n=\sum _{j=0}^{n-1}{\tau }_j(n\ge 1),$$

$$\left({\partial }_{{\tau }_n}{u}_h^{n+1},\chi \right)+A(u_h^{n+1},\chi )=(f(u_h^n),\chi )(\chi \in S_h,n=0,1,\dots ),$$

$$(w,v)$$

$$A(w,v)$$

$$⟨{\partial }_{{\tau }_n}{u}_h^{n+1},\chi ⟩+B(u_h^{n+1},\chi )=⟨f(u_h^n),\chi ⟩(\chi \in S_h,n=0,1,\dots ),$$

$$⟨w,v⟩$$

$$B(w,v)$$

$$B(w,w)=⟨w_x,w_x⟩+(2-N){\int }_I{w}_{x}wdx=|||w_x||{|}^2+\frac{N-2}{2}w(0{)}^2\ge |||w_x||{|}^2.$$

$$\mbox{$f$ is a non-decreasing function with $f(0)\geq 0$}.$$

$$\tau_{n}\geq\frac{1}{4}h^{2}{\color[rgb]{0,0,0},}$$

$$\mathcal{M}=({\mu }_{i,j}{)}_{0\le i,j\le m-1}\in {\mathbb{R}}^{m\times m},$$

$$\mathcal{A}=(a_{i,j}{)}_{0\le i,j\le m-1}\in {\mathbb{R}}^{m\times m},$$

$${\mathbf{u}}^n=(u_j^n{)}_{0\le j\le m-1}\in {\mathbb{R}}^m,$$

$${\mathbf{F}}^n=(F_j^n{)}_{0\le j\le m-1}\in {\mathbb{R}}^m,$$

$$(\mathcal{M}+{\tau }_n\mathcal{A}){\mathbf{u}}^{n+1}=\mathcal{M}{\mathbf{u}}^n+{\tau }_n{\mathbf{F}}^n(n=0,1,\dots ),$$

$$\mathcal{C}=(c_{i,j}{)}_{0\le i,j\le m-1}=\mathcal{M}+{\tau }_n\mathcal{A},c_{i,j}={\mu }_{i,j}+{\tau }_n{a}_{i,j}.$$

$${\mathbf{u}}^{n+1}={\mathcal{C}}^{-1}\left(\mathcal{M}{\mathbf{u}}^n+{\tau }_n{\mathbf{F}}^n\right)\ge 0,$$

$$\sum _{j=0}^{m-1}{c}_{i,j}>0(0\le i\le m-1).$$

$$\sum _{j=0}^{m-1}{c}_{i,j}$$

$$=\sum _{j=i-1}^{i+1}{\mu }_{i,j}+{\tau }_n{\int }_{x_{i-1}}^{x_{i+1}}{x}^{N-1}({\varphi }_{i-1}+{\varphi }_i+{\varphi }_{i+1}{)}_x({\varphi }_i{)}_{x}dx$$

$$=\sum _{j=i-1}^{i+1}{\mu }_{i,j}>0,$$

$${\tau }_n\ge -\frac{{\mu }_{i,i+1}}{a_{i,i+1}},-\frac{{\mu }_{i,i-1}}{a_{i,i-1}}(i=0,1,\cdots ,m-1),$$

Full text

Key words:

finite element method, numerical analysis, radially symmetric solution, semilinear parabolic equation

2010 Mathematics Subject Classification: 65M60, 35K58,

1 Introduction

This study was conducted to investigate the convergence property of finite element method (FEM) applied to a parabolic equation with singular coefficients for the function $u=u(x,t)$, $x\in\overline{I}=[0,1]$ , and $t\geq 0$, as expressed in

[TABLE]

where $f$ is a given locally Lipschitz continuous function, $u^{0}$ is a given continuous function, and

[TABLE]

is a given parameter.

In the study of an $N$ -dimensional semilinear heat equation, the following problem arises as

[TABLE]

where $\Omega$ represents a bounded domain in $\mathbb{R}^{N}$. If one is concerned with the radially symmetric solution $u(|\bm{x}|)=U(\bm{x})$ in the $N$ -dimensional ball \Omega=\{\bm{x}\in\mathbb{R}^{N}\mid|\bm{x}|=|\bm{x}|_{\mathbb{R}^{N}}{\color[rgb]{0,0,0}<}1\}, then (3) implies (1), where $x=|\bm{x}|$ and $u^{0}(x)=U_{0}(\bm{x})$.

For a linear case in which $f(u)=0$ is replaced by a given function $f(x,t)$, the works [7, 15] studied the convergence property of the FEM to (1) along with the corresponding steady-state problem, and two proposed schemes: the symmetric scheme, wherein they established the optimal order error estimate in the weighted $L^{2}$ norm ; and the nonsymmetric scheme, wherein they proved the $L^{\infty}$ error estimate. In this paper, both schemes are applied to the semilinear equation (1) to derive various error estimates. Moreover, this study includes a discussion of discrete positivity conservation properties, which earlier studies [7, 15] failed to embrace, but which are actually important in the study of diffusion -type equations.

Our emphasis is on FEM because we are able to use non-uniform partitions of the space variable . Therefore, the method is deemed useful for examining highly concentrated solutions at the origin. On this connection, we present our motivation for this study. The critical phenomenon appearing in the semilinear heat equation of the form

[TABLE]

in a multidimensional space has attracted considerable attention since the pioneering work of Fujita [8]. According to him, the equation is in the whole $N$ dimensional space . Any positive solution blows up in a finite time if $\alpha\leq 2/N$, whereas a solution is smooth at any time for a small initial value if $\alpha>2/N$. Therefore, expression $p_{c}=1+2/N$ is known as Fujita’s critical exponent ( [11, 6] provides some critical exponents of other equations). Generally, similar critical exponents can be found for an initial-boundary value problem for the semilinear heat equation . Some examples are given in reports of earlier studies [9, 11, 6]. However, the concrete values of those critical conditions are apparently unknown. Therefore, we found it interesting to study the numerical methods for computing the solutions of nonlinear partial differential equations in an $N$ -dimensional space. However, computing the non-stationary four-space dimensional problem is difficult, even for modern computers. We consider the FEM to solve the one space dimensional equation (1). However, we face another difficulty in dealing with the singular coefficient $(N-1)/x$, which the FEM reasonably simplified, as explained later.

As described above, the main purpose of this paper is to derive various optimal order error estimates for the symmetric and nonsymmetric schemes of [7, 15] applied to (1). These schemes are described below as (Sym) and (Non-Sym). To this end, we address mostly the general nonlinearity $f(u)$. Moreover, we study discrete positivity conservation properties. We summarize our typical results here.

• •

The solution of (Sym) is positive if $f$ and if the discretization parameters satisfy some conditions, as shown by Theorem 3.2.

• •

If $f$ is a globally Lipschitz continuous function, then the solution of (Sym) converges to the solution of (1) in the weighted $L^{2}$ norm for the space and in the $L^{\infty}$ norm for time. Moreover, the convergence is at the optimal order, as shown by Theorem 4.1.

• •

If $f$ is a locally Lipschitz continuous function and $N\leq 3$, then the solution of (Sym) converges to the solution of (1) in the weighted $L^{2}$ norm for the space and in the $L^{\infty}$ norm for time. The convergence is at the optimal order, as shown by Theorem 4.3.

• •

If $f(u)=u|u|^{\alpha}$ with $\alpha\geq 1$ and if the time partition is uniform, then the solution of (Non-Sym) converges to the solution of (1) in the $L^{\infty}(0,T;L^{\infty}(I))$ norm. The convergence is at the optimal order up to the logarithm factor, as shown by Theorem 4.6.

However, we do not proceed to applications of our schemes to the blow-up computation in this work. In fact, from the main results presented in this paper, we infer that the standard schemes of [7, 15] do not fit for the blow-up computation for large $N$. For the symmetric scheme, the restriction $N\leq 3$ reduces interest in considering radially symmetric problems. Moreover, for the nonsymmetric scheme, the use of uniform time-partitions makes it difficult to apply Nakagawa’s time-partitions control strategy: a powerful technique for computing the approximate blow-up time, as described in earlier reports [3, 12, 14, 5, 4, 13, 2]. Nevertheless, we believe that our results are of interest to researchers in this and related fields. In fact, the validity issue of the symmetric scheme only for $N\leq 3$ was pointed out earlier in [1] for a nonlinear Schrödinger equation with no mathematical evidence. The analysis reported herein reveals weak points of the two standard schemes. As a sequel to this study, we propose a new finite element scheme for (1). The scheme, which uses a nonstandard mass-lumping approximation, is shown to be positivity-preserving and convergent for any $N\geq 2$. Details will be reported in a forthcoming paper.

It is noteworthy that the finite difference method for (1) has been studied and that its optimal order convergence was proved in an earlier report [3]. Its finite difference scheme uses a special approximation around the origin to assume a uniform spatial mesh.

This paper comprises five sections. Section 2 presents our finite element schemes. Well-posedness and positivity conservation are examined in Section 3. Section 4 presents the error estimates and their proofs. Finally, Section 5 presents some numerical examples that validate our theoretical results.

2 Finite element method

First, we derive two alternate weak formulations of (1). Unless otherwise stated explicitly, we assume that $f$ is a locally Lipschitz continuous function such that

[TABLE]

Letting $\chi\in\dot{H}^{1}=\{v\in H^{1}(I)\mid v(1)=0\}$ be arbitrary, then multiplying both sides of (1a) by $x^{N-1}\chi$ and using integration by parts over $I$, we obtain

[TABLE]

Otherwise, if we multiply both sides of (1a) by $x\chi$ instead of $x^{N-1}\chi$ and integrate it over $I$, then we have

[TABLE]

We designate (4) the symmetric weak form because of the symmetric bilinear form associated with the differential operator $u_{xx}+\frac{N-1}{x}u_{x}$. In contrast, (5) is the nonsymmetric weak form. Both forms are identical at $N=2$.

We now establish the finite element schemes based on these identities. For a positive integer $m$, we introduce node points

[TABLE]

and set $I_{j}=(x_{j-1},x_{j})$ and $h_{j}=x_{j}-x_{j-1}$, where $j=1,\ldots,m$. The granularity parameter is defined as $h=\max_{1\leq j\leq m}h_{j}$. Let $\mathcal{P}_{k}(J)$ be the set of all polynomials in an interval $J$ of degree $\leq k$. We define the $\mathrm{P}1$ finite element space as

[TABLE]

Its standard basis function $\phi_{j}$, j=0,1,\cdots,m{\color[rgb]{0,0,0},} is defined as

[TABLE]

where $\delta_{ij}$ denotes Kronecker’s delta.

For time discretization, we introduce non-uniform partitions

[TABLE]

where $\tau_{j}>0$ denotes the time increments.

Generally, we write $\partial_{\tau_{n}}u_{h}^{n+1}=(u_{h}^{n+1}-u_{h}^{n})/\tau_{n}$.

We are now in a position to state the finite element schemes to be considered.

(Sym) Find $u_{h}^{n+1}\in S_{h}$, $n=0,1,\ldots$, such that

[TABLE]

where $u_{h}^{0}\in S_{h}$ is assumed to be given. Hereinafter, we set

[TABLE]

(Non-Sym) Find $u_{h}^{n+1}\in S_{h}$, $n=0,1,\ldots,$ such that

[TABLE]

where

[TABLE]

It is noteworthy that $B(\cdot,\cdot)$ is coercive in $\dot{H}^{1}$ such that

[TABLE]

3 Well-posedness and positivity conservation

In this section, we prove the following theorems.

Theorem 3.1** (Well-posedness of (Sym)).**

For a given $u_{h}^{n}\in S_{h}$ with $n\geq 0$, the scheme (Sym) admits a unique solution $u_{h}^{n+1}\in S_{h}$.

Theorem 3.2** (Positivity of (Sym)).**

In addition to the basic assumption (f1), assume that

[TABLE]

Letting $n\geq 0$ and $u_{h}^{n}\geq 0$, and assuming that

[TABLE]

then, the solution $u_{h}^{n+1}$ of (Sym) satisfies $u_{h}^{n+1}\geq 0$.

Theorem 3.3** (Comparison principle for (Sym)).**

We let $n\geq 0$ and assume that $u_{h}^{n},\tilde{u}_{h}^{n}\in S_{h}$ satisfies $u_{h}^{n}\leq\tilde{u}_{h}^{n}$ in $I$. Furthermore, we assume that (f1) and (f2) are satisfied. Similarly, we let $u_{h}^{n+1},\tilde{u}_{h}^{n+1}\in S_{h}$ be the solutions of (Sym) with $u_{h}^{n},\tilde{u}_{h}^{n}$, respectively, using the same time increment $\tau_{n}$. Moreover, we assume that (12) is satisfied. Consequently, we obtain $u_{h}^{n+1}\leq\tilde{u}_{h}^{n+1}$ in $I$ . The equality holds true if and only if $u_{h}^{n}=\tilde{u}_{h}^{n}$ in $I$.

Theorem 3.4** (Well-posedness of (Non-Sym)).**

For a given $u_{h}^{n}\in S_{h}$ with $n\geq 0$, the scheme (Non-Sym) admits a unique solution $u_{h}^{n+1}\in S_{h}$.

To prove these theorems, we conveniently rewrite (7) into a matrix form. That is, we introduce

[TABLE]

and express (7) as

[TABLE]

where $u_{m}^{n}=u_{h}^{n}(x_{m})$ is understood as $u_{m}^{n}=0$.

Lemma 3.5**.**

$\mathcal{M}$* and $\mathcal{A}$ are both tri-diagonal and positive-definite matrices.*

Theorem 3.1 is a direct consequence of this lemma. We proceed to proofs of other theorems.

Proof of Theorem 3.2.

We use the representative matrix (13) instead of (7) and set

[TABLE]

If $\mathcal{C}^{-1}\geq O$, then we obtain

[TABLE]

because $\mathcal{M}\geq O$ and $\bm{F}^{n}\geq\bm{0}$ in view of (f2). The proof that $\mathcal{C}^{-1}\geq O$ is true under (12) is divided into three steps, each described as presented below.

Step 1. We show that

[TABLE]

Letting $1\leq i\leq m-2$, we calculate

[TABLE]

because $\phi_{i-1}+\phi_{i}+\phi_{i+1}\equiv 1$ in $(x_{i-1},x_{i+1})$. Cases $i=0$ and $i=m-1$ are verified similarly.

Step 2. We show that, if

[TABLE]

then $\mathcal{C}^{-1}\geq O$. First, (15) implies that $c_{i,i-1},c_{i,i+1}\leq 0$ for $0\leq i\leq m-1$ because $a_{i,i-1},a_{i,i+1}<0$. Matrix $\mathcal{C}$ is decomposed as $\mathcal{C}=\mathcal{D}(\mathcal{I}-\mathcal{E})$, where $\mathcal{D}=(d_{i,j})_{0\leq i,j\leq m-1}$ and $\mathcal{E}=(e_{i,j})_{0\leq i,j\leq m-1}$ are defined as

[TABLE]

and where $I$ is the identity matrix. Apparently, $\mathcal{I}-\mathcal{E}$ is non-singular and $\mathcal{D}\geq O$. Using (14), we deduce

[TABLE]

Therefore, matrix $\mathcal{I}-\mathcal{E}$ is non-singular and $(\mathcal{I}-\mathcal{E})^{-1}=\sum_{k=0}^{\infty}\mathcal{E}^{k}\geq O$. Consequently, we have $\mathcal{C}^{-1}=(\mathcal{I}-\mathcal{E})^{-1}\mathcal{D}^{-1}\geq O$.

Step 3. Finally, we demonstrate that (12) implies (15). We calculate

[TABLE]

Therefore, we deduce $-\frac{\mu_{i,i+1}}{a_{i,i+1}}\leq\frac{1}{4}h^{2}$. ∎

Proof of Theorem 3.3.

Because $f(\tilde{u}_{h}^{n})-f(u_{h}^{n})\geq 0$ in $I$, the proof follows exactly the same pattern as that of the proof of Proposition 3.2. ∎

We proceed to the result for (Non-Sym):

[TABLE]

and express (9) as

[TABLE]

In view of (11), $\mathcal{M}^{\prime}$ and $\mathcal{B}$ are both tri-diagonal and positive-definite matrices. Therefore, the proof is completed.

4 Convergence and error analysis

4.1 Results

Our convergence results for (Sym) and (Non-Sym) are stated under a smoothness assumption of the solution $u$ of (1) : given $T>0$ and setting $Q_{T}=[0,1]\times[0,T]$, we assume that $u$ is sufficiently smooth such that

[TABLE]

where $\nu$ is either 0 or 1.

The partition $\{x_{i}\}_{j=0}^{m}$ of $\bar{I}=[0,1]$ is assumed to be quasi -uniform, with a positive constant $\beta$ independent of $h$ such that

[TABLE]

Finally, the approximate initial value $u_{h}^{0}$ is chosen as

[TABLE]

for a positive constant $C_{0}$.

Moreover, for $k=1,2,\ldots$, we express the positive constants $C_{k}=C_{k}(\gamma_{1},\gamma_{2},\ldots)$ and $h_{k}=h_{k}(\gamma_{1},\gamma_{2},\ldots)$ according to the parameters $\gamma_{1},\gamma_{2},\ldots$. Particularly, $C_{k}$ and $h_{k}$ are independent of $h$ and $\tau$.

Next we state the following theorems.

Theorem 4.1** (Convergence for (Sym) in $\|\cdot\|$, I).**

Assume that $f$ is a globally Lipschitz continuous function; assume (f1) and

[TABLE]

Assume that, for $T>0$, solution $u$ of (1) is sufficiently smooth that (17) for $\nu=0$ holds true. Moreover, assume that (18) and (19) are satisfied. Then, there exists $h_{1}=h_{1}(N,\beta)$ such that, for any $h\leq h_{1}$, we have

[TABLE]

where $C_{1}=C_{1}(T,M,\kappa_{0}(u),C_{0},N,\beta)$ and $u_{h}^{n}$ is the solution of (Sym).

For $L^{\infty}$ error estimates, we must further assume that $u_{h}^{0}$ is chosen as

[TABLE]

Theorem 4.2** (Convergence for (Sym) in $\|\cdot\|_{L^{\infty}(\sigma,1)}$, I).**

In addition to the assumption of Theorem 4.1, assume that (20) is satisfied. Furthermore, let $\sigma\in(0,1)$ be arbitrary. Then, there exists an $h_{2}=h_{2}(N,\beta)$ such that, for any $h\leq h_{2}$, we have

[TABLE]

where $C_{2}=C_{2}(T,M,\kappa_{0}(u),C_{0},N,\beta,\sigma)$ and $u_{h}^{n}$ is the solution of (Sym).

The restriction that $f$ is a globally Lipschitz continuous function with (f3) can be removed in the following manner.

Theorem 4.3** (Convergence of (Sym) in $\|\cdot\|$, II).**

Given that $T>0$ and that only (f1) is satisfied, we assume that (17) with $\nu=0$, (18), and (19) are satisfied. Furthermore, assume that $N\leq 3$ and that there exist positive constants $c_{1}$ and $\sigma$ such that

[TABLE]

Then there exists an $h_{3}=h_{3}(T,\kappa_{0}(u),C_{0},N,\beta)$ such that, for any $h\leq h_{3}$, we have

[TABLE]

where $C_{3}=C_{3}(T,\kappa_{0}(u),C_{0},N,\beta)$ and $u_{h}^{n}$ is the solution of (Sym).

Theorem 4.4** (Convergence for (Sym) in $\|\cdot\|_{L^{\infty}(\sigma,1)}$, II).**

Given that $T>0$ and that (f1) is satisfied, we assume that (17) with $\nu=0$, (18), (19), (20) and (21) are satisfied. Consequently, there exists $h_{4}=h_{4}(T,\kappa_{0}(u),C_{0},N,\beta)$ such that, for any $h\leq h_{4}$, we have

[TABLE]

where $C_{4}=C_{4}(T,\kappa_{0}(u),C_{0},N,\beta)$ and $u_{h}^{n}$ is the solution of (Sym).

Subsequently, let us proceed to error estimates for (Non-Sym). For the approximate initial value $u_{h}^{0}$, we choose

[TABLE]

Quasi-uniformity is also required for the time partition . Therefore, there exists a positive constant $\gamma>0$ such that

[TABLE]

where $\tau_{\min}=\min_{n\geq 0}\tau_{n}$. Moreover, we set

[TABLE]

Theorem 4.5** (Convergence for (Non-Sym), I).**

Let $f$ be a $C^{1}$ function satisfying

[TABLE]

Given $T>0$, we assume that the solution $u$ of (1) is sufficiently smooth that (17) for $\nu=1$ holds true. Furthermore, we assume that (18), (22) and (23) are satisfied. Then, there exists an $h_{5}=h_{5}(T,\kappa_{1}(u),M_{1},M_{2},\gamma,N,\beta)$ such that, for any $h\leq h_{5}$, we have

[TABLE]

where $C_{5}=C_{5}(T,\kappa_{1}(u),M_{1},M_{2},\gamma,N,\beta)>0$ and $u_{h}^{n}$ is the solution of (Non-Sym).

Finally, we state the error estimates for non-globally Lipschitz continuous function $f$. To avoid unnecessary complexity, we deal only with the power nonlinearity $f(s)=s|s|^{\alpha}$.

Theorem 4.6** (Convergence for (Non-Sym), II).**

Letting $f(s)=s|s|^{\alpha}$for $s\in\mathbb{R}$, where $\alpha\geq 1$ , then given $T>0$, we assume that (17) with $\nu=1$, (18) and (22) are satisfied. Then, there exists an $h_{6}=h_{6}(T,\kappa_{1}(u),\gamma,N,\beta)$ such that, for any $h\leq h_{6}$, we have

[TABLE]

where $C_{6}=C_{6}(T,\kappa_{1}(u),\gamma,N,\beta)$ and $u_{h}^{n}$ is the solution of (Non-Sym).

4.2 Proof of Theorems 4.1 and 4.2

We use the projection operator $P_{A}$ of $\dot{H}^{1}\to S_{h}$ associated with $A(\cdot,\cdot)$, defined for $w\in\dot{H}^{1}$ as

[TABLE]

In [7] and [10], the following error estimates are proved.

Lemma 4.7**.**

Letting $w\in C^{2}(\bar{I})\cap\dot{H}^{1}$, and (18) be satisfied, then for $h\leq h_{7}=h_{7}(N,\beta)$, we obtain

[TABLE]

where $C$ is a positive constant depending only on $N$ and $\beta$.

Proof of Theorem 4.1.

Using $P_{A}u$, we distribute the error in the form shown below.

[TABLE]

From (26), it is known that

[TABLE]

Next we derive an estimate for $\theta^{n}$. By considering the symmetric weak form (4) at $t=t_{n+1}$, we obtain

[TABLE]

which, together with (7), implies that

[TABLE]

Substituting this expression for $\chi=\theta^{n+1}$ yields the following:

[TABLE]

Correspondingly, because

[TABLE]

we provide an estimate

[TABLE]

To sum up, we obtain

[TABLE]

Therefore,

[TABLE]

where $C^{\prime}=C^{\prime}(T,\kappa_{0}(u),M,N,\beta,C_{0})>0$. By combining this expression with (28), one can deduce the desired error estimate. ∎

Proof of Theorem 4.2.

We use the same error decomposition process as that used in the previous proof where $u_{h}^{n}-u(t_{n})=\theta^{n}+\rho^{n}$ . Also, we apply (27) to estimate $\|\rho^{n}\|_{L^{\infty}(I)}$. Because

[TABLE]

we perform an estimation for $\|\theta^{n}_{x}\|$.

Substituting (29) for $\chi=\partial_{\tau_{n}}\theta^{n+1}$, we obtain the following.

[TABLE]

Correspondingly, we apply the elementary identity shown below

[TABLE]

along with Young’s inequality to obtain

[TABLE]

where $\delta_{0},\delta_{1},\delta_{2}>0$ are constants. After setting $\delta_{0}^{2}+\delta_{1}^{2}+\delta_{2}^{2}=1$, we obtain

[TABLE]

Therefore,

[TABLE]

Consequently, using (20), (30), and (31), we deduce

[TABLE]

This, together with (27) and (32), implies the desired estimate. ∎

4.3 Proof of Theorems 4.3 and 4.4

For the proof, we use the inverse inequality that follows.

Lemma 4.8** (Inverse inequality).**

Under condition (18),

[TABLE]

where $C_{\star}$ is a positive constant depending only on $N$ and $\beta$.

Proof.

Let $v_{h}\in S_{h}$ be arbitrary. From the norm equivalence in $\mathbb{R}^{2}$, we know that

[TABLE]

where $C_{\star\star}$ denotes the absolute positive constant. Given that $\|v_{h}\|_{L^{\infty}(I)}=\|v_{h}\|_{L^{\infty}(I_{1})}$, the expression is calculable as

[TABLE]

The case $\|v_{h}\|_{L^{\infty}(I)}=\|v_{h}\|_{L^{\infty}(I_{j})}$ with $j=2,\ldots,m$ is examined similarly. ∎

Proof of Theorem 4.3.

Consider (1) and (Sym) with replacement $f(s)$ in

[TABLE]

where $\mu>0$ is determined later. Then, $\tilde{f}$ satisfies condition (f3) in Theorem 4.1 such that

[TABLE]

Let $\tilde{u}$ and $\tilde{u}_{h}^{n}$ be the solutions of (1) and (Sym) with $\tilde{f}$, respectively, such that

[TABLE]

where $\theta^{n}=\tilde{u}_{h}^{n}-P_{A}\tilde{u}(t_{n})$ and $\rho^{n}=P_{A}\tilde{u}(t_{n})-\tilde{u}(t_{n})$. Applying Theorem 4.1 to $\tilde{u}$ and $\tilde{u}_{h}^{n}$, one obtains

[TABLE]

where $C_{2}=C_{2}(T,\kappa_{0}(\tilde{u}),\mu,C_{0},N,\beta)$. Moreover, an estimate (31) for $\theta^{n}$ is available. In view of Lemmas 4.7 and 4.8, we determine those estimates as

[TABLE]

where $C_{3}=C_{3}(T,\kappa_{0}(\tilde{u}),\mu,C_{0},N,\beta)$ and $C_{4}=C_{4}(N,\beta)$. Therefore, we have

[TABLE]

and

[TABLE]

At this stage, we set $\mu=1+\|u\|_{L^{\infty}(Q_{T})}$ to obtain $u=\tilde{u}$ in $Q_{T}$ by uniqueness. Moreover, because $N<4$, we can take a very small $h$ such that

[TABLE]

Consequently, $\|\tilde{u}_{h}^{n}\|_{L^{\infty}(I)}\leq\mu$ . Also, by uniqueness $u_{h}^{n}=\tilde{u}_{h}^{n}$. Therefore, (33) implies the desired conclusion. ∎

Proof of Theorem 4.4.

The proof follows the exact same pattern as that for Theorem 4.3 , but using Theorem 4.2 instead of Theorem 4.1. ∎

4.4 Proof of Theorems 4.5 and 4.6

We use the projection operator $P_{B}$ of $\dot{H}^{1}\to S_{h}$ associated with $B(\cdot,\cdot)$:

[TABLE]

In [7], the following error estimates are proved.

Lemma 4.9**.**

Letting $w\in C^{2}(\bar{I})\cap\dot{H}^{1}$ and (18) be satisfied, then for $h\leq h_{8}=h_{8}(N,\beta)$ we obtain

[TABLE]

where $C_{8}=C_{8}(N,\beta)$.

We also use a version of Poincaré’s inequality (see [15, Lemma 18.1]).

Lemma 4.10**.**

We have

[TABLE]

We can now state the proof that follows.

Proof of Theorem 4.5.

Using $P_{B}u(t)\in S_{h}$, we decompose the error into

[TABLE]

We know from (35) that

[TABLE]

Therefore, we will specifically examine estimation of $|||\theta^{n}_{x}|||$ because we are aware that

[TABLE]

Furthermore, (5) and (9) give

[TABLE]

for $\chi\in S_{h}$. Substituting this for $\chi=\theta^{n+1}$, we have

[TABLE]

This, together with (11), implies that

[TABLE]

Therefore, using (36), we deduce that

[TABLE]

These estimates actually hold . Nevertheless, their proof is postponed for Appendix A:

[TABLE]

Using (37a), (37b), (41a), and (41b), we deduce

[TABLE]

which completes the proof of Theorem 4.5. ∎

Finally, we state the following proof.

Proof of Theorem 4.6.

Consider problems (1) and (9) with replacement $f(s)=s|s|^{\alpha}$ by

[TABLE]

where $\mu>0$ is determined later. Then, $\tilde{f}$ is a $C^{1}$ function and the corresponding values of $\tilde{M}_{1}$ and $\tilde{M}_{2}$ in (f4) are expressed as $\tilde{M}_{1}=(1+\alpha)\mu^{\alpha}$ and $\tilde{M}_{2}=(1+\alpha)\alpha\mu^{\alpha-1}$.

Let $\tilde{u}$ and $\tilde{u}_{h}^{n}$ respectively represent the solutions of (1) and (9) with $\tilde{f}$. If $\mu\geq\kappa_{1}(u)$, then $u=\tilde{u}$ holds true by uniqueness. Consequently, we can apply Theorem 4.5 to obtain

[TABLE]

where $C=C(T,\kappa_{1}(u),\gamma,N,\beta)$. At this juncture, we apply small $h$ and $\tau$ such that $C\left(\log\frac{1}{h}\right)^{\frac{1}{2}}(h^{2}+\tau)<1$, and set $\mu=\kappa_{1}(u)+1$. As $\|\tilde{u}_{h}^{n}\|_{L^{\infty}(I)}\leq\kappa_{1}(u)+1=\mu$, we obtain $\tilde{u}_{h}^{n}=u_{h}^{n}$ by the uniqueness theorem. Therefore, (42) implies the desired estimate. ∎

5 Numerical examples

This section presents some numerical examples to validate our theoretical results. For this purpose, throughout this section, we set

[TABLE]

If this were the case, then the solution of (1) might blow up in the finite time. Therefore, one must devote particular attention to setting of the time increment $\tau_{n}$. Particularly, following Nakagawa [12] (see also Chen [3] and Cho–Hamada–Okamoto [5]), we use the time-increment control

[TABLE]

where $\tau=\lambda h^{2}$ and $\lambda=1/2$.

First, we compared the shapes of both solutions of (Sym) and (Non-Sym), as shown in Fig. 1 for $N=5$, $\alpha=\frac{4}{3}$ and $u(0,x)=\cos\frac{\pi}{2}x$. We used the uniform space mesh $x_{j}=jh$ ($j=0,\ldots,m$) and $h=1/m$ with $m=50$.

We computed them continuously until t_{n}{\color[rgb]{0,0,0}=}T=0.2 or $\|u_{h}\|_{2}^{-1}<\epsilon=10^{-8}$, wherein both solutions exist globally in time and approach [math] uniformly in $\overline{I}$ as $t\to\infty$. No marked differences were observed in Figs. 1(a) and  1(b). Subsequently, we took Fig. 2 for the case in which the initial value was $u(0,x)=13\cos\frac{\pi}{2}x$ . The rest of the parameters are the same. At this point, the solutions of (Sym) and (Non-Sym) blew up after $x=0.06$ with the distinct observation that the solution of the former blew up earlier than that of the latter. Furthermore, the solution of (Non-Sym) had negative values whereas that of (Sym) was always positive.

We examined the error estimates of the solutions for the same uniform space mesh $x_{j}=jh$ ($j=0,\ldots,m$) and $h=1/m$ . Also, we regarded the numerical solution with $h^{\prime}=1/480$ as the exact solution. The following quantities were compared:

[TABLE]

Fig. 3 presents results for $N=3$, $\alpha=\frac{4}{3}$ and $u(0,x)=\cos\frac{\pi}{2}x$. We used the uniform time increment $\tau_{n}=\tau=\lambda h^{2}$ $(n=0,1,\ldots)$ with $\lambda=1/2$ and computed until $t\leq T=0.005$. For (Sym), we observed the theoretical convergence rate $h^{2}+\tau$ in the $\|\cdot\|$ norm (see Theorem 4.3) , whereas the rate in the $L^{\infty}$ norm deteriorated slightly. For (Non-Sym), we observed second-order convergence in the $L^{\infty}$ norm, which supports the results presented in Theorem 4.5.

Moreover, we considered the case for $N=4$, which is not supported in Theorem 4.3 for (Sym) . Also, we chose $\alpha=4$ and $u(0,x)=3\cos\frac{\pi}{2}x$ for this case. Fig. 4(d) displays the shape of the solution, which blew up at approximately $T=0.0035$. Furthermore, we computed errors until $T=0.0011,0.0022$, and $0.0033$ using the uniform meshes $x_{j}$ and $\tau_{n}$ with $\lambda=0.11$. From Fig. 4, we observed the second-order convergence in the $\|\cdot\|$ norm, suggesting the possibility of removing assumption $N\leq 3$.

Finally, we observed the non-increasing property of the energy functional. The energy functional associated with (1) is given as

[TABLE]

We can use the standard method to prove that $J(t)$ is non-increasing in $t$.

This non-increasing property plays an important role in the blow-up analysis of the solution of (1) , as presented by Nakagawa [12]. Therefore, it is of interest whether a discrete version of this non-increasing property holds true. Actually, introducing the discrete energy functional associated with (Sym) as

[TABLE]

we prove the following. Appendix B presents the proof.

Proposition 5.1**.**

$J_{h}(n)$* is a non-increase sequence of $n$.*

Now let $N=3$, $\alpha=\frac{4}{3}$, and $u(0,x)=\cos\frac{\pi}{2}x,~{}13\cos\frac{\pi}{2}x$. We determined the time increment $\tau_{n}$ through (43) for the uniform space mesh $x_{j}=jh$ with $h=1/m$ and $m=50$. Fig. 5 presents the results, which support that of Proposition 5.1.

Appendix A Proofs of (41a) and (41b)

Proofs of (41a) and (41b) are stated in this appendix using the same notation as that used in Section 4.

Proof of (41a).

By application of (39), (37a), and (37b), we derived the expression

[TABLE]

Consequently, we have

[TABLE]

Therefore, similarly to the derivation of (31), we obtain from (22) the expression of

[TABLE]

to complete the proof. ∎

Proof of (41b).

First, we prove the case of $n=0$. Substituting (38) for $n=0$ and $\chi=\theta^{1}$, we obtain

[TABLE]

Because $\theta^{0}=0$, we apply (37b) to get

[TABLE]

Repeatedly using $\theta^{0}=0$, we obtain

[TABLE]

Next we assume $n\geq 0$ and $t_{n+2}\leq T$. Consequently, from (38), we derive

[TABLE]

for any $\chi\in S_{h}$. Substituting this expression for $\chi=\partial_{\tau_{n+1}}\theta^{n+2}$, we obtain

[TABLE]

Here, we accept the following estimates:

[TABLE]

In view of the quasi-uniformity of time partition (23), we have

[TABLE]

Summing up, we deduce

[TABLE]

where $b_{n}=|||\partial_{\tau_{n}}\theta^{n+1}|||$. Therefore,

[TABLE]

which, together with (44), implies the desired inequality (41b).

We now prove (46a)–(46c).

*Estimation for $J_{1}$. * We apply Taylor’s theorem to obtain

[TABLE]

where $s_{1}=u_{h}^{n+1}-\mu_{1}(u_{h}^{n+1}-u_{h}^{n})$ and $s_{2}=u(t_{n+1})-\mu_{2}[u(t_{n+1})-u(t_{n})]$ for some $\mu_{1},\mu_{2}\in[0,1]$. In view of (37a), (37b), and (41a), we find the following estimates

[TABLE]

and

[TABLE]

*Estimation for $J_{2}$. * We begin with

[TABLE]

where $s_{3}=u(t_{n+1})+\mu_{3}(u(t_{n+2})-u(t_{n+1}))$ and $s_{4}=u(t_{n+1})+\mu_{4}(u(t_{n})-u(t_{n+1}))$ for some $\mu_{3},\mu_{4}\in[0,1]$, $\eta_{1}\in[t_{n+1},t_{n+2}]$, and $\eta_{2}\in[t_{n},t_{n+1}]$. Next, we obtain the following estimate :

[TABLE]

[TABLE]

*Estimation for $J_{3}$. * We express $J_{3}$ as

[TABLE]

for some $s_{5}\in[t_{n+1},t_{n+2}]$, $s_{6}\in[t_{n},t_{n+1}]$ and $s_{7}\in[s_{6},s_{5}]\subset[t_{n},t_{n+2}]$. Therefore,

[TABLE]

*Estimation for $J_{4}$. * For some $s_{8}\in[t_{n+1},t_{n+2}]$, $s_{9}\in[t_{n},t_{n+1}]$ , and $s_{10}\in[s_{9},s_{8}]$, we obtain the expression

[TABLE]

Therefore, using (35),

[TABLE]

∎

Appendix B Proof of Proposition 5.1

Proof.

Substituting $\chi=\partial_{\tau_{n}}u_{h}^{n+1}$ for (7), we have

[TABLE]

Therefore, for the conditions

[TABLE]

we obtain

[TABLE]

which implies that $J_{h}(n+1)\leq J_{h}(n)$.

We can validate (48a) and (48b) . Also, (48a) is derived readily. To prove (48b), we set $g(s)=\frac{1}{\alpha+2}|s|^{\alpha+2}$, and apply the mean value theorem to deduce

[TABLE]

where $w=w(x)=u_{h}^{n}+\sigma(u_{h}^{n+1}-u_{h}^{n})$ and $\sigma=\sigma(x)\in(0,1)$. Consequently,

[TABLE]

Then we repeat the mean value theorem to resolve

[TABLE]

where $\tilde{w}=u^{n}_{h}+\tilde{\sigma}(w-u_{h}^{n})$ and $\tilde{\sigma}=\tilde{\sigma}(x)\in(0,1)$. Therefore,

[TABLE]

which gives (48b). ∎

Acknowledgments.

This work was supported by JST CREST Grant No. JPMJCR15D1, Japan, and JSPS KAKENHI Grant No. 15H03635, Japan. In addition, the first author was supported by the Program for Leading Graduate Schools, MEXT, Japan.