Optimal error estimates of a linearized second-order BDF scheme for a nonlocal parabolic problem

Abstract

Purpose

This paper focuses on the unconditionally optimal error estimates of a linearized second-order scheme for a nonlocal nonlinear parabolic problem. The first step of the scheme is based on Crank–Nicholson method while the second step is the second-order BDF method.

Design/methodology/approach

A rigorous error analysis is done, and optimal L² error estimates are derived using the error splitting technique. Some numerical simulations are presented to confirm the study’s theoretical analysis.

Findings

Optimal L² error estimates and energy norm.

Originality/value

The goal of this research article is to present and establish the unconditionally optimal error estimates of a linearized second-order BDF finite element scheme for the reaction-diffusion problem. An optimal error estimate for the proposed methods is derived by using the temporal-spatial error splitting techniques, which split the error between the exact solution and the numerical solution into two parts, that is, the temporal error and the spatial error. Since the spatial error is not dependent on the time step, the boundedness of the numerical solution in L∞-norm follows an inverse inequality immediately without any restriction on the grid mesh.

Citation

Daoussa Haggar, M.S. and Mbehou, M. (2024), "Optimal error estimates of a linearized second-order BDF scheme for a nonlocal parabolic problem", Arab Journal of Mathematical Sciences, Vol. 30 No. 1, pp. 112-129. https://doi.org/10.1108/AJMS-05-2022-0126

1. Introduction

In this paper, we consider the following parabolic problem with nonlocal nonlinearity:

(1.1)$\left\{\begin{array}{c}{u}_t-a\left(l\left(u\right)\right)\mathrm{\Delta }u+\alpha |u{|}^{p-2}u=f\left(u\right)in\mathrm{\Omega }\times \left(0,T\right],\hfill \\ u\left(x,t\right)=0\text{on}\partial \mathrm{\Omega }\times \left(0,T\right],\hfill \\ u\left(x,0\right)=u_0\left(x\right)\text{in}\mathrm{\Omega },\hfill \end{array}\right.$

where $\mathrm{\Omega }\subset {\mathbb{R}}^d$, d ≥ 1 is again a domain with a smooth boundary ∂Ω, a and f are functions to be defined in the next section and l denote a continuous linear form on L²(Ω) given by

$l\left(u\left(t\right)\right)={\int }_{\mathrm{\Omega }}g\left(x\right)u\left(t,x\right)dx,$

where g is a function on L²(Ω).

The study of nonlocal parabolic problems has received considerable attention in recent years ([1–3] and the references therein). This kind of problems arises in various situations, for instance, u could describe the density of a population (for instance, bacteria) subject to spreading. The diffusion coefficient a is then supposed to depend on the entire population in the domain rather than on the local density, that is, moves are guided by considering the global state of the medium. The problem is nonlocal in the sense that the diffusion coefficient is determined by a global quantity. Besides its mathematical motivation because of the presence of the nonlocal term a(l(u)), such problems come from physical situations related to migration of a population of bacteria in a container in which the velocity of migration v = a∇u depends on the global population in a subdomain Ω′ ⊂ Ω given by a(l(u)).

Simsen and Ferreira [4] have discussed not only the existence and uniqueness of solutions for this problem but also continuity with respect to initial values, the exponential stability of weak solutions and important results on the existence of a global attractor. The numerical methods for the nonlocal problems have been investigated by many authors as like in Refs [5, 6] and the references therein. However, they are restricted to nonlocal reaction terms or nonlocal boundary conditions. Chaudhary et al. [7] studied the convergence analysis of the Crank–Nicolson finite element method for the nonlocal problem involving the Dirichlet energy. Mbehou et al. [8] studied (1.1) using the Crank–Nicolson Galerkin finite element method. The main focus on this paper was to present the exponential decay and vanishing of the solutions in finite time. They also derived the optimal convergence order in L²-norm using P_r with r ≥ 1 finite elements. Yin and Xu [9] applied the finite-volume method to obtain approximate solutions for a nonlocal problem on reactive flows in porous media and derived the optimal convergence order in the L² norm. Almeida et al. [10] presented convergence analysis for a fully discretized approximation to a nonlocal problem involving a parabolic equation with moving boundaries, with the finite element method applied for the space variables and the Crank–Nicolson method for the time. Recently, Yang et al. [11] derived the unconditional optimal error estimate of Galerkin FEMs for the time-dependent Klein–Gordon–Schrodinger equations using the error splitting technique. Also in Ref. [12], Yang and Jiang applied the linearized second-order backward differentiation formulae (BDF) Galerkin Finite element methods (FEMs) for the Landau-Lifshitz equations to derive the unconditional optimal error estimates.

Our goal in this research article is to give and establish the unconditionally optimal error estimates of a linearized second-order BDF finite element scheme for the reaction-diffusion problem (1.1). Using P_r (r ≥ 1) finite element to approximate the solution of (1.1), the optimal error estimates O(Δt² + h^r+1) in L² norm are derived using the error splitting technique.

This paper is organized as follows. In Section 2, we recall few known results and present few regularities, which are used in the proof of the optimal error estimates. To prove the optimal error estimates by the error splitting technique, the temporal errors and the spatial errors are shown in Sections 3 and 4, respectively. Finally numerical results are presented in Section 5 to demonstrate our theoretical analysis.

2. Preliminaries and main results

Let $\mathrm{\Omega }\subset {\mathbb{R}}^d$ (d ≥ 1) be a bounded domain with a smooth boundary ∂Ω = Γ. The standard notations (see for instance Refs [13, 14]) will be used throughout this work. The Lebesgue space is denoted L^p(Ω), 1 ≤ p ≤ ∞, with norms $\Vert \cdot {\Vert }_{L^p}$ but the L²(Ω)-norm will be denoted by ‖ ⋅‖. For any nonnegative integer m and real number p ≥ 1, the classical Sobolev spaces:

$W^{m,p}\left(\mathrm{\Omega }\right)=\left\{v\in L^p\left(\mathrm{\Omega }\right);D^{\alpha }v\in L^p\left(\mathrm{\Omega }\right)\text{for\hspace{0.17em}all}\left|\alpha \right|\le m\right\},$

equipped with the semi-norm

$|v{|}_{m,p}={\left\{\sum _{\left|\alpha \right|=m}{\int }_{\mathrm{\Omega }}|D^{\alpha }v{|}^{p}dx\right\}}^{1/p},$

and the norm

(2.1)$\Vert v{\Vert }_{m,p}={\left\{\sum _{0\le |\alpha |\le m}{\int }_{\mathrm{\Omega }}|D^{\alpha }v{|}^{p}dx\right\}}^{1/p},$

with the usual extension when p = ∞. When p = 2, W^m,p(Ω) is the Hilbert space H^m(Ω) with the scalar product:

${\left(\left(v,w\right)\right)}_m=\sum _{|\alpha |\le m}\left(D^{\alpha }v,D^{\alpha }w\right)$

The norm of H^m(Ω) will be denoted by ‖ ⋅‖_m. It should be mentioned that D^α stands for the derivative in the sense of distribution, while α = (α₁, …, α_d) denotes a multi-index of length |α| = α₁ + ⋯ + α_d. We also employ the standard notation of Bochner spaces, such as L^q(0, T, X) with norm

$\begin{array}{l}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

where X is an Hilbert space and ‖ ⋅‖_X the norm of X. For all these notions on Sobolev spaces and Bochner spaces, we refer to Refs [13, 15].

Throughout this paper, the following known inequalities will be frequently used [13].

(2.2)$\Vert v{\Vert }_{L^r}\le C\Vert v{\Vert }_1\left(2\le r\le 6\right)\forall v\in H^1\left(\mathrm{\Omega }\right)$

(2.3)$\Vert v{\Vert }_{\infty }\le C\Vert v{\Vert }^{1/2}\Vert v{\Vert }_2^{1/2}\forall v\in H^2\left(\mathrm{\Omega }\right)\text{.}$

Let us now suppose that α is a nonnegative constant and p > 1. Simsen and Ferreira [4] proved the existence and uniqueness of global solution under the following hypotheses.

H1.

u₀ ∈ L²(Ω).

H2.

$f:\mathbb{R}\to \mathbb{R}$ is Lipschitz–continuous function, that is, there exists γ > 0 such that |f(s) − f(t)| ≤ γ|s − t|, for all$s,t\in \mathbb{R}$ and f(0) = 0.

H3.

$a:\mathbb{R}\to \mathbb{R}$ is bounded with 0 < m ≤ a(s) ≤ M, for all $s\in \mathbb{R}$ with ${\lambda }_1>\frac{\gamma }{m}$, where λ₁ is the first eigenvalue of $\left(-\mathrm{\Delta },H_0^1\left(\mathrm{\Omega }\right)\right)$.

H4.

$a:\mathbb{R}\to \mathbb{R}$ is Lipschitz–continuous with $|a\left(s_1\right)-a\left(s_2\right)|\le A|s_1-s_2|,\forall s_1,s_2\in \mathbb{R}$.

H5.

$l:L^2\left(\mathrm{\Omega }\right)\to \mathbb{R}$ is a continuous linear form, i.e. there exists g ∈ L²(Ω) such that l(u) = l_g(u) = ∫_Ω g(x)u(x)dx, for all u ∈ L²(Ω).

Theorem 2.1

(Existence and uniqueness of solution, [4]). Assume that p ≥ 2 and if the hypotheses (H1)–(H5) hold, then problem (1.1) possesses a unique solution, that is, there exists a unique function u such that

(2.4) $u\in L^2\left(0,T,H_0^1\left(\mathrm{\Omega }\right)\cap L^p\left(\mathrm{\Omega }\right)\right)\cap \mathcal{C}\left(\left[0,T\right];L^2\left(\mathrm{\Omega }\right)\right),$

(2.5)$u_t\hspace{0.17em}\in L^2\left(0,T,H^{-1}\left(\mathrm{\Omega }\right)\right),$

(2.6) $u\left(x,0\right)=u_0\left(x\right)in\mathrm{\Omega }.$

(2.7) $\left(u_t,v\right)+a\left(l\left(u\right)\right)\left(\mathrm{\nabla }u,\mathrm{\nabla }v\right)+\alpha \left(|u{|}^{p-2}u,v\right)=\left(f\left(u\right),v\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right),$

where (2.7)must be understood as an equality in ${\mathcal{D}}^{\prime }\left(0,T\right)$.

Given the hypotheses (H1)–(H5), we will also adopt another hypothesis, namely

H6.

for all r ≥ 1,

(2.8) $\Vert u_0{\Vert }_{H^{r+1}}+\Vert u{\Vert }_{L^{\infty }\left(H^{r+1}\left(\mathrm{\Omega }\right)\right)}+\Vert u_t{\Vert }_{L^2\left(H^{r+1}\left(\mathrm{\Omega }\right)\right)}+\Vert u_{tt}{\Vert }_{L^2\left(H^1\left(\mathrm{\Omega }\right)\right)}+\Vert u_{\mathit{ttt}}{\Vert }_{L^2\left(L^2\left(\mathrm{\Omega }\right)\right)}\le C.$

The following lemmas will be useful.

Lemma 2.1

(cf. Ref. [16]). For all p ∈ (1, ∞) and τ ≥ 0, there exists a generic constant C = C(p, d) such that for all $\xi ,\eta \in {\mathbb{R}}^d$ with d ≥ 1 we have

(2.9) $\Vert \xi {|}^{p-2}\xi -|\eta {|}^{p-2}\eta |\le C|\xi -\eta {|}^{1-\tau }{\left(|\xi |+|\eta |\right)}^{p-2+\tau }.$

(2.10)$\left(|\xi {|}^{p-2}\xi -|\eta {|}^{p-2}\eta \right)\cdot \left(\xi -\eta \right)\ge C|\xi -\eta {|}^{2+\tau }{\left(|\xi |+|\eta |\right)}^{p-2-\tau }.$

Lemma 2.2

(cf. Ref. [3]). Let a and b be two nonnegaitve numbers. Then for all s ∈ (1, ∞),

(2.11) $|a^s-b^s|\le |a-b|{\left(a+b\right)}^{s-1}.$

Lemma 2.3

(cf. Ref. [17]). Let a_k, b_k, c_k and γ_k, for integers k ≥ 0, be the positive numbers such that

(2.12) $a_n+\tau \sum _{k=0}^n{b}_k\le \tau \sum _{k=0}^n{\gamma }_k{a}_k+\tau \sum _{k=0}^n{c}_k+B,\text{for}n\ge 0.$

Suppose that τγ_k < 1, for all k, and set ${\sigma }_k=\frac{1}{\left(1-\tau {\gamma }_k\right)}$. Then

(2.13)$a_n+\tau \sum _{k=0}^n{b}_k\le \exp \left(\tau \sum _{k=0}^n{\gamma }_k{\sigma }_k\right)\left(\tau \sum _{k=0}^n{c}_k+B\right),\text{for}n\ge 0.$

Remark. If the first sum on the right hand side of (2.12) extends only up to n − 1, then estimate (2.13) holds for all k > 0 with σ_k = 1.

Lemma 2.4

(H^k-estimate of elliptic equations [18]). Suppose that v is a solution of the boundary value problem

$\begin{array}{l}\hfill -\mathrm{\Delta }v=f,\text{in}\mathrm{\Omega },\\ \hfill v=0,\text{on}\partial \mathrm{\Omega },\end{array}$

where $\mathrm{\Omega }\subset {\mathbb{R}}^d$, d = 2, 3, is a smooth and bounded domain. Then,

(2.14)$\Vert v{\Vert }_{H^k}\le C\Vert f{\Vert }_{H^{k-2}},k=2,3.$

Let ${\mathcal{T}}_h=\left\{K\right\}$ be a uniform triangular or tetrahedral partition of Ω into triangles or tetrahedrons. Thus, let $h={\max }_{K\in {\mathcal{T}}_h}\left\{h_K\right\}$ denote the mesh size, where h_K = diam(K) = max{‖x − y‖, x, y ∈ K}, and V_h be the finite dimensional subspace of $H_0^1\left(\mathrm{\Omega }\right)$, which consists of continuous piecewise polynomials of degree r ≥ 1 on ${\mathcal{T}}_h$.

Let {t_n| t_n = nΔt; 0 ≤ n ≤ N} be a uniform partition of [0, T] with time step Δt = T/N. We write uⁿ = u(x, t_n), Uⁿ ≈ u(x, t_n) and for any sequence of functions ${\left\{w^n\right\}}_{n=0}^N$ define

$\begin{array}{ll}\hfill D_1{w}^n& =\frac{w^n-w^{n-1}}{\mathrm{\Delta }t},\hfill \\ \hfill {\overline{w}}^n& =\frac{1}{2}\left(w^n+w^{n-1}\right),n=\mathrm{1,2}\dots ,N,\hfill \\ \hfill D_2{w}^n& =\frac{3}{2}{D}_1{w}^n-\frac{1}{2}{D}_1{w}^{n-1}\text{and}{\widehat{w}}^n=2w^{n-1}-w^{n-2},n=2,\dots ,N.\hfill \end{array}$

The following telescope formula is for n ≥ 2

(2.15)$\left(D_2{w}^n,w^n\right)=\frac{1}{4\mathrm{\Delta }t}\left(\Vert w^n{\Vert }^2-\Vert w^{n-1}{\Vert }^2+\Vert {\widehat{w}}^{n+1}{\Vert }^2-\Vert {\widehat{w}}^n{\Vert }^2+\Vert w^n-2w^{n-1}+w^{n-2}{\Vert }^2\right).$

Under the above notations, we propose the following linearized second-order BDF Galerkin finite element scheme associated to (1.1), which is to find $U_h^n\in V_h$ such that

• Step 1: For $U_h^0={\mathrm{\Pi }}_h{u}_0\in V_h$, find $U_h^1\in V_h$ such that for all v_h ∈ V_h

(2.16)$(D_1{U}_h^1,v_h)+a(l({\widehat{U}}_h^1))(\mathrm{\nabla }{\overline{U}}_h^1,\mathrm{\nabla }{v}_h)+\alpha (|{\widehat{U}}_h^1{|}^{p-2}{\overline{U}}_h^1,v_h)=(f({\widehat{U}}_h^1),v_h),$

where ${\widehat{U}}_h^1$ is given by

(2.17)$\left(\frac{{\widehat{U}}_h^1-U_h^0}{\mathrm{\Delta }t/2},v_h\right)+a(l\left(U_h^0\right)(\mathrm{\nabla }{\widehat{U}}_h^1,\mathrm{\nabla }{v}_h)+\alpha (|U_h^0{|}^{p-2}{\widehat{U}}_h^1,v_h)=\left(f\left(U_h^0\right),v_h\right)\text{.}$

• Step 2: For 2 ≤ n ≤ N, find $U_h^n\in V_h$ such that for all v_h ∈ V_h

(2.18)$(D_2{U}_h^n,v_h)+a(l({\widehat{U}}_h^n))(\mathrm{\nabla }{U}_h^n,\mathrm{\nabla }{v}_h)+\alpha (|{\widehat{U}}_h^n{|}^{p-2}{U}_h^n,v_h)=(f({\widehat{U}}_h^n),v_h)\text{.}$

Π_h is an interpolation operator from $H_0^1\left(\mathrm{\Omega }\right)$ to V_h.

Theorem 2.2

Assume that the hypotheses (H1)– (H5) hold. Then the fully discrete system defined in (2.16)–(2.18)has a unique solution $U_h^n$ which satisfies

(2.19)$\Vert U_h^n{\Vert }^2+\Vert {\widehat{U}}_h^{n+1}{\Vert }^2+C\mathrm{\Delta }t\sum _{k=1}^n\Vert \mathrm{\nabla }{U}_h^k{\Vert }^2\le C\Vert u_0{\Vert }^2.$

Proof. 1 For the existence, taking $v_h=U_h^1$, $v_h={\widehat{U}}_h^1$ and $v_h=U_h^n$ in (2.16)–(2.18), respectively, the existence and uniqueness of $U_h^1$, ${\widehat{U}}_h^1$ and $U_h^n$ are from the Lax–Milgram theorem and the hypothesis (H3).

Let $v_h={\overline{U}}_h^n$ in (2.16), we have

$\frac{1}{2\mathrm{\Delta }t}\left(\Vert U_h^1{\Vert }^2-\Vert U_h^0{\Vert }^2\right)+a(l({\widehat{U}}_h^1))\Vert \mathrm{\nabla }{\overline{U}}_h^1{\Vert }^2+\alpha (|{\widehat{U}}_h^1{|}^{p-2}{\overline{U}}_h^1,{\overline{U}}_h^1)=(f({\widehat{U}}_h^1),{\overline{U}}_h^1)\le \Vert f({\widehat{U}}_h^1)\Vert \Vert {\overline{U}}_h^1\Vert .$

Drop the third term of the left hand side, use the lower bound of a(⋅) and (H2),

$\frac{1}{2\mathrm{\Delta }t}\left(\Vert U_h^1{\Vert }^2-\Vert U_h^0{\Vert }^2\right)+M\Vert \mathrm{\nabla }{\overline{U}}_h^1{\Vert }^2\le L\Vert {\widehat{U}}_h^1\Vert \Vert {\overline{U}}_h^1\Vert \le C\Vert {\widehat{U}}_h^1{\Vert }^2+\frac{M}{2}\Vert \mathrm{\nabla }{\overline{U}}_h^1{\Vert }^2.$

and

(2.20)$\Vert U_h^1{\Vert }^2+C\mathrm{\Delta }t\Vert \mathrm{\nabla }{\overline{U}}_h^1{\Vert }^2\le C\mathrm{\Delta }t\Vert {\widehat{U}}_h^1{\Vert }^2+\Vert u_0{\Vert }^2.$

Now, let $v_h={\widehat{U}}_h^1$ in (2.17), the same arguments used above give us

$\Vert {\widehat{U}}_h^1{\Vert }^2+C\mathrm{\Delta }t\Vert \mathrm{\nabla }{\widehat{U}}_h^1{\Vert }^2\le C\Vert u_0{\Vert }^2.$

Taking $v_h=U_h^n$ in (2.18), using the lower bound of a(⋅), (H2) and dropping the third term of the left hand side lead to

$\left(D_2{U}_h^n,U_h^n\right)+M\Vert \mathrm{\nabla }{U}_h^n{\Vert }^2\le L\Vert {\widehat{U}}_h^n\Vert \Vert U_h^n\Vert .$

From the telescope (2.15), we obtain

$\begin{array}{l}\hfill \frac{1}{4\mathrm{\Delta }t}(\Vert U_h^n{\Vert }^2-\Vert U_h^{n-1}{\Vert }^2+\Vert {\widehat{U}}_h^{n+1}{\Vert }^2-\Vert {\widehat{U}}_h^n{\Vert }^2)+M\Vert \mathrm{\nabla }{U}_h^n{\Vert }^2\le CL\Vert {\widehat{U}}_h^n\Vert \Vert \mathrm{\nabla }{U}_h^n\Vert \\ \hfill \le C\Vert {\widehat{U}}_h^n{\Vert }^2+\frac{M}{2}\Vert \mathrm{\nabla }{U}_h^n{\Vert }^2.\end{array}$

That is

(2.21)$\Vert U_h^n{\Vert }^2-\Vert U_h^{n-1}{\Vert }^2+\Vert {\widehat{U}}_h^{n+1}{\Vert }^2-\Vert {\widehat{U}}_h^n{\Vert }^2+\frac{M}{2}\mathrm{\Delta }t\Vert \mathrm{\nabla }{U}_h^n{\Vert }^2\le C\Vert {\widehat{U}}_h^n{\Vert }^2.$

The relation (2.19) is obtained by summing up the above relation (2.21) and using the discrete Gronwall lemma 2.3. 

The main result of this work is presented in the following theorem.

Theorem 2.3

Suppose that system (1.1)has a unique solution u satisfying (H6). Then the fully discrete system defined in (2.16)–(2.18)has a unique solution $U_h^n$, and

(2.22) $\underset{0\le n\le N}{\max }(\Vert u^n-U_h^n{\Vert }^2+\Vert {\widehat{u}}^n-{\widehat{U}}_h^n{\Vert }^2+\mathrm{\Delta }t\sum _{k=0}^n\Vert \mathrm{\nabla }\left(u^n-U_h^n\right){\Vert }^2)\le C\left(\mathrm{\Delta }{t}^4+h^{2r+2}\right),$

where C is a positive constant independent of Δt and h.

The proof of this theorem will be done in the following sections.

3. Error estimates for the semi-discrete problem

Let us introduce the corresponding time discrete system associated with (1.1)

• Step 1: for U⁰ = u₀, find U¹ by

(3.1)$\left\{\begin{array}{c}{D}_1{U}^1-a(l({\widehat{U}}^1))\mathrm{\Delta }{\overline{U}}^1+\alpha |{\widehat{U}}^1{|}^{p-2}{\overline{U}}^1=f({\widehat{U}}^1)\hfill \\ U^1=0\text{on}\partial \mathrm{\Omega },\hfill \end{array}\right.$

where ${\widehat{U}}^1$ is the solution to

(3.2)$\left\{\begin{array}{c}\frac{{\widehat{U}}^1-U^0}{\mathrm{\Delta }t/2}-a(l(U^0))\mathrm{\Delta }{\widehat{U}}^1+\alpha |U^0{|}^{p-2}{\widehat{U}}^1=f(U^0)\hfill \\ {\widehat{U}}^1=0\text{on}\partial \mathrm{\Omega }.\hfill \end{array}\right.$

• Step 2: for 2 ≤ n ≤ N, find Uⁿ by

(3.3)$\left\{\begin{array}{c}{D}_2{U}^n-a(l({\widehat{U}}^n))\mathrm{\Delta }{U}^n+\alpha |{\widehat{U}}^n{|}^{p-2}{U}^n=f({\widehat{U}}^n)\hfill \\ U^n=0\text{on}\partial \mathrm{\Omega }.\hfill \end{array}\right.$

The weak formulations of (3.1)–(3.3) are defined as follows: find $U^n\in H_0^1(\mathrm{\Omega })$ such that for all $v\in H_0^1(\mathrm{\Omega })$

(3.4)$(D_1{U}^1,v)+a(l({\widehat{U}}^1))(\mathrm{\nabla }{\overline{U}}^1,\mathrm{\nabla }v)+\alpha (|{\widehat{U}}^1{|}^{p-2}{\overline{U}}^1,v)=(f({\widehat{U}}^1),v),$

and for 2 ≤ n ≤ N

(3.5)$(D_2{U}^n,v)+a(l({\widehat{U}}^n))(\mathrm{\nabla }{U}^n,\mathrm{\nabla }v)+\alpha (|{\widehat{U}}^n{|}^{p-2}{U}^n,v)=(f({\widehat{U}}^n),v),$

with ${\widehat{U}}^1\in H_0^1(\mathrm{\Omega })$ such that

(3.6)$\left(\frac{{\widehat{U}}^1-U^0}{\mathrm{\Delta }t/2},v\right)+a(l(U^0))(\mathrm{\nabla }{\widehat{U}}^1,\mathrm{\nabla }v)+\alpha (|U^0{|}^{p-2}{\widehat{U}}^1,v)=(f(U^0),v)\text{.}$

The existence and uniqueness of the solution to problems (3.4)–(3.6) can be easily proved by using Lax–Milgram theorem.

Let u be the exact solution of (1.1). Then, u satisfies the following equations:

(3.7)$D_1{u}^1-a(l({\widehat{u}}^1))\mathrm{\Delta }{\overline{u}}^1+\alpha |{\widehat{u}}^1{|}^{p-2}{\overline{u}}^1=f({\widehat{u}}^1)+R^1$

(3.8)$D_2{u}^n-a(l({\widehat{u}}^n))\mathrm{\Delta }{u}^n+\alpha |{\widehat{u}}^n{|}^{p-2}{u}^n=f({\widehat{u}}^n)+R^n,n=2,\dots ,N,$

where ${\widehat{u}}^1$ satisfies

(3.9)$\frac{{\widehat{u}}^1-u_0}{\mathrm{\Delta }t/2}-a\left(l\left(u_0\right)\right)\mathrm{\Delta }{\widehat{u}}^1+\alpha |u_0{|}^{p-2}{\widehat{u}}^1=f\left(u_0\right)+{\widehat{R}}^0.$

${\widehat{R}}^0$, R¹ and Rⁿ are, respectively, the truncation errors given by

${\widehat{R}}^0=\left(\frac{{\widehat{u}}^1-u_0}{\mathrm{\Delta }t/2}-u_t^{1/2}\right)-\mathrm{\Delta }{\widehat{u}}^1(a(l(u_0))-a(l(u^{1/2})))+\alpha |u_0{|}^{p-2}(u^{1/2}-u_0)+\alpha (|u_0{|}^{p-2}{u}_0-|u^{1/2}{|}^{p-2}{u}^{1/2})+(f(u_0)-f(u^{1/2})),R^1=(D_{1}1u^1-u_t^{1/2})-\mathrm{\Delta }{\overline{u}}^1(a(l({\widehat{u}}^1))-a(l(u^{1/2})))-a(l(u^{1/2}))\mathrm{\Delta }({\overline{u}}^1-u^{1/2})+\alpha |{\widehat{u}}^1{|}^{p-2}({\overline{u}}^1-{\widehat{u}}^1)+\alpha (|{\widehat{u}}^1{|}^{p-2}{\widehat{u}}^1-|u^{1/2}{|}^{p-2}{u}^{1/2})+(f({\widehat{u}}^1)-f(u^{1/2}))R^n=(D_2{u}^n-u_t^n)-\mathrm{\Delta }{u}^n(a(l({\widehat{u}}^n))-a(l(u^n)))+\alpha |{\widehat{u}}^n{|}^{p-2}(u^n-{\widehat{u}}^n)+\alpha (|{\widehat{u}}^n{|}^{p-2}{\widehat{u}}^n-|u^n{|}^{p-2}{u}^n)+(f({\widehat{u}}^n)-f(u^n))\text{.}$

By Taylor formula and relation (2.9) with τ = 1, it is easy to see that

(3.10)${\left(\sum _{n=1}^N\mathrm{\Delta }t\Vert R^n{\Vert }_1^2\right)}^{1/2}+\mathrm{\Delta }t\Vert {\widehat{R}}^0{\Vert }_1\le C\mathrm{\Delta }{t}^2.$

Let us denote

${\overline{e}}^1={\overline{u}}^1-{\overline{U}}^1;{\widehat{e}}^n={\widehat{u}}^n-{\widehat{U}}^n;e^n=u^n-U^n\text{for}1\le n\le N.$

We have the following assumption.

Lemma 3.1

Assume that the exact solution u of (1.1)satisfies the regularities (2.8). Then there exists a positive constant C independent of Δt such that

(3.11)$\Vert {\widehat{e}}^1\Vert +\mathrm{\Delta }{t}^{1/2}\Vert \mathrm{\nabla }{\widehat{e}}^1\Vert \le C\mathrm{\Delta }{t}^2.$

Proof. Subtracting (3.6) from (3.9) leads to

$2{\widehat{e}}^1-\mathrm{\Delta }ta\left(l\left(u_0\right)\right)\mathrm{\Delta }{\widehat{e}}^1+\alpha \mathrm{\Delta }t|u_0{|}^{p-2}{\widehat{e}}^1=\mathrm{\Delta }t{\widehat{R}}^0.$

Testing the above equation by ${\widehat{e}}^1$ yield

$2\Vert {\widehat{e}}^1{\Vert }^2+\mathrm{\Delta }ta\left(l\left(u_0\right)\right)\Vert \mathrm{\nabla }{\widehat{e}}^1{\Vert }^2+\alpha \mathrm{\Delta }t{\int }_{\mathrm{\Omega }}|u_0{|}^{p-2}|{\widehat{e}}^1{|}^{2}dx=\mathrm{\Delta }t({\widehat{R}}^0,{\widehat{e}}^1)$

Using the left bound of a(⋅) to the left hand side and Young’s inequality to the right hand side, we obtain

$2\Vert {\widehat{e}}^1{\Vert }^2+m\mathrm{\Delta }t\Vert \mathrm{\nabla }{\widehat{e}}^1{\Vert }^2+\alpha \mathrm{\Delta }t{\int }_{\mathrm{\Omega }}|u_0{|}^{p-2}|{\widehat{e}}^1{|}^{2}dx\le \frac{1}{2}\mathrm{\Delta }{t}^2\Vert {\widehat{R}}^0{\Vert }^2+\frac{1}{2}\Vert {\widehat{e}}^1{\Vert }^2.$

The proof ended by dropping the third term of the left hand side and applying (3.10) to the right hand side.

Based upon (3.11), we have

Proposition 3.1

Suppose that the solution u of (1.1)satisfies the regularities (2.8). Then there exists a generic constant C that does not dependent on Δt such that

(3.12) $\Vert e^1{\Vert }^2+\mathrm{\Delta }t\Vert \mathrm{\nabla }{e}^1{\Vert }^2\le C\mathrm{\Delta }{t}^5.$

Proof. Subtracting (2.16) from (3.7) and observing that e⁰ = 0 leads to

$\begin{array}{l}\hfill e^1-\mathrm{\Delta }ta(l(u_0))\mathrm{\Delta }{\overline{e}}^1+\alpha \mathrm{\Delta }t|{\widehat{u}}^1{|}^{p-2}{\overline{e}}^1=-\alpha \mathrm{\Delta }t{\overline{u}}^1(|{\widehat{u}}^1{|}^{p-2}-|{\widehat{U}}^1{|}^{p-2})\\ \hfill +\mathrm{\Delta }t\mathrm{\Delta }{\overline{u}}^1(a(l({\widehat{u}}^1))-a(l({\widehat{U}}^1)))+\mathrm{\Delta }t(f({\widehat{u}}^1)-f({\widehat{U}}^1))+\mathrm{\Delta }t{R}^1.\end{array}$

Testing the above equation by e¹ and using the fact that ${\overline{e}}^1=\frac{1}{2}{e}^1$, we have

(3.13)$\begin{array}{l}\Vert e^1{\Vert }^2+\frac{1}{2}m\mathrm{\Delta }t\Vert \mathrm{\nabla }{e}^1{\Vert }^2+\frac{1}{2}\alpha \mathrm{\Delta }t{\int }_{\mathrm{\Omega }}|{\widehat{u}}^1{|}^{p-2}|e^1{|}^{2}dx\le \alpha \mathrm{\Delta }t({\overline{u}}^1(|{\widehat{u}}^1{|}^{p-2}-|{\widehat{U}}^1{|}^{p-2}),e^1)\\ +\mathrm{\Delta }t(a(l({\widehat{u}}^1))-a(l({\widehat{U}}^1)))(\mathrm{\nabla }{\overline{u}}^1,\mathrm{\nabla }{e}^1)+\mathrm{\Delta }t(f({\widehat{u}}^1)-f({\widehat{U}}^1),e^1)+\mathrm{\Delta }t(R^1,e^1)\hfill \\ =\sum _{i=1}^4{I}_i.\hfill \end{array}$

We have

$\begin{array}{ll}{I}_1& =\alpha \mathrm{\Delta }t({\overline{u}}^1(|{\widehat{u}}^1{|}^{p-2}-|{\widehat{U}}^1{|}^{p-2}),e^1)\le C_1(\Vert {\overline{u}}^1{\Vert }_{L^{\infty }},\Vert {\widehat{u}}^1{\Vert }_{L^{\infty }},\Vert {\widehat{U}}^1{\Vert }_{L^{\infty }},p)\mathrm{\Delta }t\Vert {\widehat{e}}^1\Vert \Vert e^1\Vert \hfill \\ \text{mtd}\end{array}$

$\begin{array}{ll}{I}_2& =\mathrm{\Delta }t(a(l({\widehat{u}}^1))-a(l({\widehat{U}}^1)))(\mathrm{\nabla }{\overline{u}}^1,\mathrm{\nabla }{e}^1)\le \mathrm{\Delta }tA\Vert {\widehat{e}}^1\Vert \Vert \mathrm{\nabla }{\overline{u}}^1\Vert \Vert \mathrm{\nabla }{e}^1\Vert (\text{using\hspace{0.17em}}(\mathbf{H}\mathbf{4}))\hfill \\ \text{mtd}\end{array}$

$\begin{array}{ll}{I}_3& =\mathrm{\Delta }t(f({\widehat{u}}^1)-f({\widehat{U}}^1),e^1)\le \gamma \mathrm{\Delta }t\Vert {\widehat{e}}^1\Vert \Vert e^1\Vert (\text{using\hspace{0.17em}}(\mathbf{H}\mathbf{2}))\hfill \\ \text{mtd}\end{array}$

$I_4=\mathrm{\Delta }t(R^1,e^1)\le C_4\Vert R^1{\Vert }^2+\frac{1}{4}\Vert e^1{\Vert }^2\le C\mathrm{\Delta }{t}^5+\frac{1}{4}\Vert e^1{\Vert }^2(\text{using\hspace{0.17em}}(3.11))\text{.}$

Taking these estimates into (3.13), we obtain the desire result.

The main result in this section is as follows.

Theorem 3.1

Suppose that the solution u of (1.1)satisfies the regularities (2.8). Then there exists a generic constant C that does not dependent on Δt such that

(3.14) $\underset{1\le n\le N}{\max }\left(\Vert e^n{\Vert }^2+\Vert {\widehat{e}}^n{\Vert }^2+\mathrm{\Delta }t\sum _{k=1}^n\Vert \mathrm{\nabla }{e}^k{\Vert }^2\right)\le C\mathrm{\Delta }{t}^4,$

(3.15)$\underset{0\le n\le N}{\max }\Vert U^n{\Vert }_{\infty }\le C,$

where C is a positive constant independent of n and Δt.

Proof. The proof of this theorem will be done using the mathematical induction.

In view of (3.11) and (3.12), the inequality (3.14) holds for n = 0, 1. Since U⁰ = u₀, the inquality (3.15) holds for n = 0. Now, let us assume that (3.14) and (3.15) hold for n ≤ m with m ≤ N − 1. Then we need to prove the inequality for n = m + 1. By the definition of ${\widehat{U}}^n$ and the induction assumption, $\Vert {\widehat{U}}^n{\Vert }_{\infty }\le C$.

Subtracting (2.18) from (3.8) results in the following equation:

(3.16)$\begin{array}{l}{D}_2{e}^n-a(l({\widehat{U}}^n))\mathrm{\Delta }{e}^n+\alpha |{\widehat{U}}^n{|}^{p-2}{e}^n=(a(l({\widehat{u}}^n))-a(l({\widehat{U}}^n)))\mathrm{\Delta }{u}^n\\ -\alpha u^n(|{\widehat{u}}^n{|}^{p-2}-|{\widehat{U}}^n{|}^{p-2})+(f({\widehat{u}}^n)-f({\widehat{U}}^n))+R^n.\hfill \end{array}$

Multiply (3.17) by 4Δteⁿ and integrate it over Ω. The use of the telescope formula to the resulting equation leads to

$\begin{array}{ll}\Vert e^n{\Vert }^2& -\Vert e^{n-1}{\Vert }^2+\Vert {\widehat{e}}^{n+1}{\Vert }^2-\Vert {\widehat{e}}^n{\Vert }^2+\Vert e^n-2e^{n-1}+e^{n-2}{\Vert }^2+4\mathrm{\Delta }ta(l({\widehat{U}}^n))\Vert \mathrm{\nabla }{e}^n{\Vert }^2\hfill \\ \text{mtd}\\ \text{mtd}\end{array}$

Use the lower bound of a(⋅) and drop certain positive terms on the left hand side of the above equation leads to

(3.17)$\begin{array}{ll}\Vert e^n{\Vert }^2& -\Vert e^{n-1}{\Vert }^2+\Vert {\widehat{e}}^{n+1}{\Vert }^2-\Vert {\widehat{e}}^n{\Vert }^2+4m\mathrm{\Delta }t\Vert \mathrm{\nabla }{e}^n{\Vert }^2\le 4\mathrm{\Delta }t(a(l({\widehat{u}}^n))-a(l({\widehat{U}}^n)))(\mathrm{\nabla }{u}^n,\mathrm{\nabla }{e}^n)\hfill \\ \text{mtd}\\ \text{mtd}\end{array}$

We have

$\begin{array}{ll}{J}_1=4\mathrm{\Delta }t(a(l({\widehat{u}}^n))-a(l({\widehat{U}}^n)))(\mathrm{\nabla }{u}^n,\mathrm{\nabla }{e}^n)& \le 4\mathrm{\Delta }tA\Vert {\widehat{e}}^n\Vert \Vert \mathrm{\nabla }{u}^n\Vert \Vert \mathrm{\nabla }{e}^n\Vert (\text{using\hspace{0.17em}}(\mathbf{H}\mathbf{4}))\hfill \\ \text{mtd}\end{array}$

$\begin{array}{ll}{J}_2=4\mathrm{\Delta }t\alpha (u^n(|{\widehat{u}}^n{|}^{p-2}-|{\widehat{U}}^n{|}^{p-2}),e^n)& \le C(\Vert u^n{\Vert }_{L^{\infty }},\Vert {\widehat{u}}^n{\Vert }_{L^{\infty }},\Vert {\widehat{U}}^n{\Vert }_{L^{\infty }},p)\mathrm{\Delta }t\Vert {\widehat{e}}^n\Vert \Vert e^n\Vert \hfill \\ \text{mtd}\end{array}$

$\begin{array}{ll}{J}_3& =4\mathrm{\Delta }t(f({\widehat{u}}^n)-f({\widehat{U}}^n),e^n)\le 4\gamma \mathrm{\Delta }t\Vert {\widehat{e}}^n\Vert \Vert e^n\Vert (\text{using\hspace{0.17em}}(\mathbf{H}\mathbf{2}))\hfill \\ \text{mtd}\end{array}$

$J_4=4\mathrm{\Delta }t(R^n,e^n)\le C\mathrm{\Delta }t\left(\Vert R^n{\Vert }^2+\Vert e^n{\Vert }^2\right).$

Therefore,

$\Vert e^n{\Vert }^2-\Vert e^{n-1}{\Vert }^2+\Vert {\widehat{e}}^{n+1}{\Vert }^2-\Vert {\widehat{e}}^n{\Vert }^2+\mathrm{\Delta }t\Vert \mathrm{\nabla }{e}^n{\Vert }^2\le C\mathrm{\Delta }t\left(\Vert e^n{\Vert }^2+\Vert {\widehat{e}}^n{\Vert }^2+\Vert R^n{\Vert }^2\right).$

Summing up the above inequality and using the discrete Gronwall inequality, we get

$\Vert e^n{\Vert }^2+\Vert {\widehat{e}}^n{\Vert }^2+\mathrm{\Delta }t\sum _{k=1}^n\Vert \mathrm{\nabla }{e}^k{\Vert }^2\le C\mathrm{\Delta }{t}^4.$

From ‖eⁿ‖ ≤ CΔt², we have

$\begin{array}{l}\text{mtd}\\ \text{mtd}\end{array}$

Applying Lemma 2.4 for the linear elliptic problems (3.1)–(3.3) with the induction assumptions gives the H² estimate

$\begin{array}{ll}\Vert U^n{\Vert }_2& \le C\Vert D_k{U}^n\Vert +C\Vert \mathrm{\nabla }{U}^n\Vert +C\Vert |{\widehat{U}}^n{|}^{p-2}{U}^n\Vert +L\Vert {\widehat{U}}^n\Vert \hfill \\ \text{mtd}\end{array}$

Using (2.3), we have

$\Vert U^n{\Vert }_{\infty }\le C\Vert U^n{\Vert }^{1/2}\Vert U^n{\Vert }_2^{1/2}\le C$

which concludes the proof.

4. Error estimates for the fully discrete problem

In this section, we will prove the optimal spatial error estimates. Let Π_h be an interpolation operator and $R_h:H_0^1(\mathrm{\Omega })\to V_h$ be a Ritz projection operator defined by

(4.1)${\int }_{\mathrm{\Omega }}\mathrm{\nabla }(u-R_{h}u)\cdot \mathrm{\nabla }wdx=0,\forall w\in H_0^1(\mathrm{\Omega })\text{.}$

Then we have the following lemma.

Lemma 4.1

(cf. Ref. [19]). If $u\in H^{r+1}(\mathrm{\Omega })\cap H_0^1(\mathrm{\Omega })$, then

(4.2) $\Vert u-{\mathrm{\Pi }}_{h}u\Vert +h\Vert \mathrm{\nabla }(u-{\mathrm{\Pi }}_{h}u)\Vert \le C{h}^{r+1}\Vert u{\Vert }_{H^{r+1}},$

(4.3) $\Vert u-R_{h}u\Vert +h\Vert \mathrm{\nabla }(u-R_{h}u)\Vert \le C{h}^{r+1}\Vert u{\Vert }_{H^{r+1}}$

(4.4) $\Vert v_h{\Vert }_{\infty }\le C{h}^{-d/2}\Vert v_h\Vert ,\forall v_h\in V_h.$

where C is a positive constant that does not depend on h and r.

Let us denote

$\begin{array}{l}\text{mtd}\\ \text{mtd}\\ \text{mtd}\end{array}$

From lemma 4.1, we have

(4.5)$\Vert E_h^0\Vert +h\Vert \mathrm{\nabla }{E}_h^0\Vert \le C{h}^{r+1}$

(4.6)$\Vert {\overline{E}}^1\Vert +\Vert {\widehat{E}}^n\Vert +\Vert E^n\Vert +h(\Vert \mathrm{\nabla }{\overline{E}}^1\Vert +\Vert \mathrm{\nabla }{\widehat{E}}^n\Vert +\Vert \mathrm{\nabla }{E}^n\Vert )\le C{h}^{r+1}.$

Lemma 4.2

Assume that the exact solution u of (1.1)satisfies the regularities (2.8). Then there exists a positive constant C independent of Δt and h such that

(4.7) $\Vert {\widehat{e}}_h^1{\Vert }^2+\mathrm{\Delta }t\Vert \mathrm{\nabla }{\widehat{e}}_h^1{\Vert }^2\le C{h}^{2r+2}$

(4.8)$\Vert {\widehat{U}}_h^1\Vert \le C.$

Proof. From equations (2.17) and (3.6), ${\widehat{e}}_h^1$ satisfies the following equation:

$\begin{array}{l}2({\widehat{e}}_h^1,v_h)+\mathrm{\Delta }ta(l\left(\left(U_h^0\right)\right)(\mathrm{\nabla }{\widehat{e}}_h^1,\mathrm{\nabla }{v}_h)+\alpha \mathrm{\Delta }t(|U_h^0{|}^{p-2}{\widehat{e}}_h^1,v_h)=2\left(E_h^0,v_h\right)-2({\widehat{E}}^1,v_h)\\ -\mathrm{\Delta }t(a(l\left(\left(U^0\right)\right)-a\left(l\left(\left(U_h^0\right)\right)\right)(\mathrm{\nabla }{\widehat{U}}^1,\mathrm{\nabla }{v}_h)-\alpha \mathrm{\Delta }t(|U_h^0{|}^{p-2}{\widehat{E}}^1,v_h)\\ -\alpha \mathrm{\Delta }t(\left(|U^0{|}^{p-2}-|U_h^0{|}^{p-2}\right){\widehat{U}}^1,v_h)+\mathrm{\Delta }t\left(f\left(U^0\right)-f\left(U_h^0\right),v_h\right)\text{.}\end{array}$

Setting $v_h={\widehat{e}}_h^1$ in the above equations leads to

(4.9)$\begin{array}{ll}2\Vert {\widehat{e}}_h^1{\Vert }^2& +m\mathrm{\Delta }t\Vert \mathrm{\nabla }{\widehat{e}}_h^1{\Vert }^2+\alpha \mathrm{\Delta }t{\int }_{\mathrm{\Omega }}|U_h^0{|}^{p-2}|{\widehat{e}}_h^1{|}^{2}dx\le 2(E_h^0+{\widehat{E}}^1,{\widehat{e}}_h^1)\hfill \\ \text{mtd}\\ \text{mtd}\\ \text{mtd}\end{array}$

From (4.5) and (4.6), we have

$\begin{array}{ll}{L}_1& =2(E_h^0+{\widehat{E}}^1,{\widehat{e}}_h^1)\le C{h}^{2r+2}+\frac{1}{2}\Vert {\widehat{e}}_h^1{\Vert }^2,\hfill \\ L_2& =\mathrm{\Delta }t(a(l\left(\left(U^0\right)\right)-a\left(l\left(\left(U_h^0\right)\right)\right)(\mathrm{\nabla }{\widehat{U}}^1,\mathrm{\nabla }{\widehat{e}}_h^1)\le \mathrm{\Delta }tA\Vert E_h^0\Vert \Vert \mathrm{\nabla }{\widehat{U}}^1\Vert \Vert \mathrm{\nabla }{\widehat{e}}_h^1\Vert \hfill \\ \text{mtd}\\ L_3& =\alpha \mathrm{\Delta }t(|U_h^0{|}^{p-2}{\widehat{E}}^1,{\widehat{e}}_h^1)\le \alpha \mathrm{\Delta }tC\left(\Vert U_h^0{\Vert }_{\infty }\right)\Vert {\widehat{E}}^1\Vert \Vert {\widehat{e}}_h^1\Vert \le C\mathrm{\Delta }{t}^2{h}^{2r+2}+\frac{1}{4}\Vert {\widehat{e}}_h^1{\Vert }^2,\hfill \\ L_4& =\alpha \mathrm{\Delta }t(\left(|U^0{|}^{p-2}-|U_h^0{|}^{p-2}\right){\widehat{U}}^1,{\widehat{e}}_h^1)\le C(\Vert U^0{\Vert }_{\infty },\Vert U_h^0{\Vert }_{\infty },\Vert {\widehat{U}}^1{\Vert }_{\infty })\mathrm{\Delta }t\Vert E_h^0\Vert \Vert {\widehat{e}}_h^1\Vert \hfill \\ \text{mtd}\\ L_5& =\mathrm{\Delta }t(f\left(U^0\right)-f\left(U_h^0\right),{\widehat{e}}_h^1)\le \gamma \mathrm{\Delta }t\Vert E_h^0\Vert \Vert {\widehat{e}}_h^1\Vert \le C\mathrm{\Delta }{t}^2{h}^{2r+2}+\frac{1}{4}\Vert {\widehat{e}}_h^1{\Vert }^2.\hfill \end{array}$

Combining these estimates into (4.9), we get (4.7).

From the inverse inequality, $\Vert {\widehat{e}}_h^1{\Vert }_{\infty }\le C{h}^{-d/2}\Vert {\widehat{e}}_h^1\Vert \le C{h}^{r/2}$ and

$\Vert {\widehat{U}}_h^1{\Vert }_{\infty }\le \Vert R_h{U}^1{\Vert }_{\infty }+\Vert {\widehat{e}}_h^1{\Vert }_{\infty }\le C.$

The main result in this section is as follows.

Theorem 4.1

Suppose that the exact solution u of (1.1)satisfies the regularities (2.8). Then there exists a positive constant C independent of Δt and h such that

(4.10) $\underset{1\le n\le N}{\max }\left(\Vert e_h^n{\Vert }^2+\Vert {\widehat{e}}_h^n{\Vert }^2+\mathrm{\Delta }t\sum _{k=1}^n\Vert \mathrm{\nabla }{e}_h^k{\Vert }^2\right)\le C{h}^{2r+2},$

(4.11) $\underset{0\le n\le N}{\max }\Vert U_h^n\Vert \le C.$

Proof. The proof of this result will be done by mathematical induction. Since $U_h^0={\mathrm{\Pi }}_h{U}^0$, (4.11) holds for n = 0. To compute the error estimate (4.10) for n = 1, subtract (3.4) from (2.16) and take $v_h={\overline{e}}_h^1\equiv \frac{1}{2}(e_h^1+E_h^0)$,

(4.12)$\begin{array}{ll}\Vert e_h^1{\Vert }^2& +2m\mathrm{\Delta }t\Vert \mathrm{\nabla }{\overline{e}}_h^1{\Vert }^2+2\alpha \mathrm{\Delta }t{\int }_{\mathrm{\Omega }}|{\widehat{U}}_h^1{|}^{p-2}|{\overline{e}}_h^1{|}^{2}dx\le (E^1+E_h^0,{\overline{e}}_h^1)\hfill \\ \text{mtd}\\ \text{mtd}\\ \text{mtd}\end{array}$

From (4.5) and (4.6), we have

$\begin{array}{ll}{T}_1& =(E^1+E_h^0,{\overline{e}}_h^1)\le C{h}^{2r+2}+\frac{1}{8}\Vert e_h^1{\Vert }^2,\hfill \\ T_2& =2\mathrm{\Delta }t(a(l(({\widehat{U}}^1))-a(l(({\widehat{U}}_h^1)))(\mathrm{\nabla }{\overline{U}}^1,\mathrm{\nabla }{\overline{e}}_h^1)\le \mathrm{\Delta }tA\Vert {\widehat{U}}^1-{\widehat{U}}_h^1\Vert \Vert \mathrm{\nabla }{\overline{U}}^1\Vert \Vert \mathrm{\nabla }{\overline{e}}_h^1\Vert \hfill \\ \text{mtd}\\ T_3& =2\alpha \mathrm{\Delta }t(|{\widehat{U}}_h^1{|}^{p-2}{\overline{E}}^1,{\overline{e}}_h^1)\le \alpha \mathrm{\Delta }tC(\Vert {\widehat{U}}_h^1{\Vert }_{\infty })\Vert {\overline{E}}^1\Vert \Vert {\overline{e}}_h^1\Vert \le C\mathrm{\Delta }{t}^2{h}^{2r+2}+\frac{1}{8}\Vert e_h^1{\Vert }^2,\hfill \\ T_4& =2\alpha \mathrm{\Delta }t((|{\widehat{U}}^1{|}^{p-2}-|{\widehat{U}}_h^1{|}^{p-2}){\overline{U}}^1,{\overline{e}}_h^1)\le C(\Vert {\widehat{U}}^1{\Vert }_{\infty },\Vert {\widehat{U}}_h^1{\Vert }_{\infty },\Vert {\overline{U}}^1{\Vert }_{\infty })\mathrm{\Delta }t\Vert {\widehat{U}}^1-{\widehat{U}}_h^1\Vert \Vert {\overline{e}}_h^1\Vert \hfill \\ \text{mtd}\\ T_5& =2\mathrm{\Delta }t(f({\widehat{U}}^1)-f({\widehat{U}}_h^1),{\overline{e}}_h^1)\le \gamma \mathrm{\Delta }t\Vert {\widehat{U}}^1-{\widehat{U}}_h^1\Vert \Vert {\overline{e}}_h^1\Vert \le C\mathrm{\Delta }{t}^2{h}^{2r+2}+C\mathrm{\Delta }{t}^2\Vert {\widehat{e}}_h^1{\Vert }^2+\frac{1}{8}\Vert e_h^1{\Vert }^2.\hfill \end{array}$

Taking these estimates into (4.12) and using Lemma 4.2, we conclude that

(4.13)$\Vert e_h^1{\Vert }^2+\mathrm{\Delta }t\Vert \mathrm{\nabla }{\overline{e}}_h^1{\Vert }^2\le C{h}^{2r+2}$

which proves (4.10) for n = 1.

Now, we assume that (4.10) and (4.11) hold for n = m − 1, 2 ≤ m ≤ N, then we need to show it also holds for n = m. By the definition of ${\widehat{U}}_h^n$ and the induction assumption, $\Vert {\widehat{U}}_h^n{\Vert }_{\infty }\le C$.

Subtracting (2.18) from (3.5), we obtain

$\begin{array}{l}(D_2{e}_h^n,v_h)+a(l(({\widehat{e}}_h^n))(\mathrm{\nabla }{e}_h^n,\mathrm{\nabla }{v}_h)+\alpha (|{\widehat{e}}_h^n{|}^{p-2}{e}_h^n,v_h)=-(D_2{E}^n,v_h)\\ -(a(l(({\widehat{U}}^n))-a(l(({\widehat{U}}_h^n)))(\mathrm{\nabla }{U}^n,\mathrm{\nabla }{v}_h)-\alpha (|{\widehat{U}}_h^n{|}^{p-2}{E}^n,v_h)\\ -\alpha ((|{\widehat{U}}^n{|}^{p-2}-|{\widehat{U}}_h^n{|}^{p-2})U^n,v_h)+(f({\widehat{U}}^n)-f({\widehat{U}}_h^n),v_h)\text{.}\end{array}$

If one takes $v_h=4\mathrm{\Delta }t{e}_h^n$ and uses the telescope formula, one obtains

$\begin{array}{ll}\Vert e_h^n{\Vert }^2& -\Vert e_h^{n-1}{\Vert }^2+\Vert {\widehat{e}}_h^{n+1}{\Vert }^2-\Vert {\widehat{e}}_h^n{\Vert }^2+\Vert e_h^n-2e_h^{n-1}+e_h^{n-2}{\Vert }^2+4\mathrm{\Delta }ta(l({\widehat{U}}_h^n))\Vert \mathrm{\nabla }{e}_h^n{\Vert }^2\hfill \\ \text{mtd}\\ \text{mtd}\end{array}$

That is

(4.14)$\begin{array}{ll}\Vert e_h^n{\Vert }^2& -\Vert e_h^{n-1}{\Vert }^2+\Vert {\widehat{e}}_h^{n+1}{\Vert }^2-\Vert {\widehat{e}}_h^n{\Vert }^2+4m\mathrm{\Delta }t\Vert \mathrm{\nabla }{e}_h^n{\Vert }^2\le 4\mathrm{\Delta }t(D_2{E}^n,e_h^n)\hfill \\ \text{mtd}\\ \text{mtd}\\ \text{mtd}\end{array}$

The quantities K_i, i = 1, …, 5 can be bounded by the similar way T_i, i = 1, …, 5:

$\begin{array}{ll}{K}_1& =4\mathrm{\Delta }t(D_2{E}^n,e_h^n)\le C\mathrm{\Delta }t(\Vert D_2{E}^n{\Vert }^2+\Vert e_h^n{\Vert }^2)\le C\mathrm{\Delta }t{h}^{2r+2}+C\mathrm{\Delta }t\Vert e_h^n{\Vert }^2,\hfill \\ K_2& =4\mathrm{\Delta }t(a(l(({\widehat{U}}^n))-a(l(({\widehat{U}}_h^n)))(\mathrm{\nabla }{U}^n,\mathrm{\nabla }{e}_h^n)\le \mathrm{\Delta }tA\Vert {\widehat{E}}^n+{\widehat{e}}_h^n\Vert \Vert \mathrm{\nabla }{U}^n\Vert \Vert \mathrm{\nabla }{e}_h^n\Vert \hfill \\ \text{mtd}\\ K_3& =4\alpha \mathrm{\Delta }t(|{\widehat{U}}_h^n{|}^{p-2}{E}^n,e_h^n)\le 4\alpha \mathrm{\Delta }tC(\Vert {\widehat{U}}_h^n{\Vert }_{\infty })\Vert E^n\Vert \Vert e_h^n\Vert \le C\mathrm{\Delta }t{h}^{2r+2}+C\mathrm{\Delta }t\Vert e_h^n{\Vert }^2,\hfill \\ K_4& =4\alpha \mathrm{\Delta }t((|{\widehat{U}}^n{|}^{p-2}-|{\widehat{U}}_h^n{|}^{p-2})U^n,e_h^n)\le C(\Vert {\widehat{U}}^n{\Vert }_{\infty },\Vert {\widehat{U}}_h^n{\Vert }_{\infty },\Vert U^n{\Vert }_{\infty })\mathrm{\Delta }t\Vert {\widehat{E}}^n+{\widehat{e}}_h^n\Vert \Vert e_h^n\Vert \hfill \\ \text{mtd}\\ K_5& =4\mathrm{\Delta }t(f({\widehat{U}}^n)-f({\widehat{U}}_h^n),e_h^n)\le 4\gamma \mathrm{\Delta }t\Vert {\widehat{E}}^n+{\widehat{e}}_h^n\Vert \Vert e_h^n\Vert \le C\mathrm{\Delta }t{h}^{2r+2}+C\mathrm{\Delta }t(\Vert {\widehat{e}}_h^n{\Vert }^2+\Vert e_h^n{\Vert }^2)\text{.}\hfill \end{array}$

Taking these bounds into (4.14), we obtain

$\Vert e_h^n{\Vert }^2-\Vert e_h^{n-1}{\Vert }^2+\Vert {\widehat{e}}_h^{n+1}{\Vert }^2-\Vert {\widehat{e}}_h^n{\Vert }^2+4m\mathrm{\Delta }t\Vert \mathrm{\nabla }{e}_h^n{\Vert }^2\le C\mathrm{\Delta }t{h}^{2r+2}+C\mathrm{\Delta }t(\Vert {\widehat{e}}_h^n{\Vert }^2+\Vert e_h^n{\Vert }^2)\text{.}$

Sum up the above inequality and use the discrete Gronwall Lemma 2.3 leads to

$\Vert e_h^n{\Vert }^2+\Vert {\widehat{e}}_h^n{\Vert }^2+\mathrm{\Delta }t\sum _{k=1}^n\Vert \mathrm{\nabla }{e}_h^k{\Vert }^2\le \exp (CT)h^{2r+2}.$

From the inverse inequality, $\Vert e_h^n{\Vert }_{\infty }\le C{h}^{-d/2}\Vert e_h^n\Vert \le C{h}^{r/2}$ and

$\Vert U_h^n{\Vert }_{\infty }\le \Vert R_h{U}^n{\Vert }_{\infty }+\Vert e_h^n{\Vert }_{\infty }\le C.$

5. Numerical results

In this section, we present several numerical simulations to illustrate our theoretical analysis. Since the resulting matrix of the linear system (2.16)–(2.18) is sparse, symmetric and positive definite, an incomplete Cholesky factorization is performed and the result is used as preconditioner in the preconditioned conjugate method iterative solver (see for instance Refs [20, 21]).

To analyze the convergence rate, we consider the following problem.

(5.1)$\left\{\begin{array}{c}{u}_t-a(l(u))\mathrm{\Delta }u+\alpha |u{|}^{p-2}u=f(u)+g\text{in}\mathrm{\Omega }\times (0,T]\hfill \\ u(x,t)=0\text{on}\partial \mathrm{\Omega }\times (0,T]\hfill \\ u(x,0)=u_0(x)\text{in}\mathrm{\Omega },\hfill \end{array}\right.$

with Ω = (0,1)², the coefficient α = 1, p = 3.5

$f(s)=s(10-s),a(s)=3+\cos (s),l(u)={\int }_{\mathrm{\Omega }}udx.$

g is chosen correspondingly to the exact solution

$u(x,y,t)=2(1+t^{2}exp(-t))xy(1-x)(1-y).$

We simulated the above problem on uniform meshes with a linear finite element approximation (r = 1) and T = 0.1.

For the convergence with respect to the mesh size h, we choose Δt = h² and we solve problem (2.16)(2.18) with different values of h (h = 1/5; 1/10; 1/15; 1/20; 1/25); from our theoretical analysis, the L²-norm errors are in order O(h² + Δt²) = O(h² + h⁴) ∼ O(h²). H¹-norm errors are in order O(h + Δt²) = O(h + h⁴) ∼ O(h). In Figure 1, we plot the log of errors against log(h). One can see that for L²-norm, the slope is almost 2, and for H¹ − norm, the slope is almost 1, which are in good agreement with our theoretical analysis.

For the convergence with respect to the time step Δt, h is fixed (h = 0.01), and we solve problem (2.16)(2.18) with different time steps Δt = 0.1; 0.05; 0.025; 0.0125 (Δt = 0.1 × 2^1−l, l = 1, …, 4), and the L²-norm errors are in order O(h² + Δt²) ∼ O(Δt²). Figure 2 shows the plots of log  L²-error norm against log(Δt). Again, one can see that the slope is almost 2. These results are consistent with our theoretical analysis.

6. Conclusion

We have presented and analyzed a linearized second-order BDF Galerkin finite element method for the nonlocal parabolic problems. We have proved the L² and energy error estimates using sufficient conditions on the exact solution. We also presented some numerical experiments on Matlab’s environment, and our numerical results confirm the theoretical analysis. The results in this paper lay the foundation for developing finite element based numerical methods for more general and complicated nonlocal problems both stationary and evolutionary.

Figures

Figure 1

Convergence rate with respect to the mesh size h in L² and H¹ norm

Figure 2

Convergence rate with respect to the time step Δt in L² norm