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)ut−a(l(u))Δu+α|u|p−2u=f(u)inΩ×(0,T],u(x,t)=0on∂Ω×(0,T],u(x,0)=u0(x)inΩ,

where Ω⊂Rd, 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(u(t))=∫Ωg(x)u(t,x)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 Ω⊂Rd (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 ‖⋅‖Lp but the L²(Ω)-norm will be denoted by ‖ ⋅‖. For any nonnegative integer m and real number p ≥ 1, the classical Sobolev spaces:

Wm,p(Ω)={v∈Lp(Ω);Dαv∈Lp(Ω)for all|α|≤m},

equipped with the semi-norm

|v|m,p=∑|α|=m∫Ω|Dαv|pdx1/p,

and the norm

(2.1)‖v‖m,p=∑0≤|α|≤m∫Ω|Dαv|pdx1/p,

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

((v,w))m=∑|α|≤m(Dαv,Dαw)

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

‖w‖Lq(X)=∫0T‖w(t)‖Xqdt1/q,1≤q<∞,‖w‖L∞(X)=esssup0≤t≤T‖w(t)‖X,

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)‖v‖Lr≤C‖v‖1(2≤r≤6)∀v∈H1(Ω)

(2.3)‖v‖∞≤C‖v‖1/2‖v‖21/2∀v∈H2(Ω).

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:R→R is Lipschitz–continuous function, that is, there exists γ > 0 such that |f(s) − f(t)| ≤ γ|s − t|, for alls,t∈R and f(0) = 0.

H3.

a:R→R is bounded with 0 < m ≤ a(s) ≤ M, for all s∈R with λ1>γm, where λ₁ is the first eigenvalue of (−Δ,H01(Ω)).

H4.

a:R→R is Lipschitz–continuous with |a(s1)−a(s2)|≤A|s1−s2|,∀s1,s2∈R.

H5.

l:L2(Ω)→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∈L2(0,T,H01(Ω)∩Lp(Ω))∩C([0,T];L2(Ω)),

(2.5)ut ∈L2(0,T,H−1(Ω)),

(2.6) u(x,0)=u0(x)inΩ.

(2.7) (ut,v)+a(l(u))(∇u,∇v)+α(|u|p−2u,v)=(f(u),v),∀v∈H01(Ω),

where (2.7) must be understood as an equality in D′(0,T).

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

H6.

for all r ≥ 1,

(2.8) ‖u0‖Hr+1+‖u‖L∞(Hr+1(Ω))+‖ut‖L2(Hr+1(Ω))+‖utt‖L2(H1(Ω))+‖uttt‖L2(L2(Ω))≤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 ξ,η∈Rd with d ≥ 1 we have

(2.9) ‖ξ|p−2ξ−|η|p−2η|≤C|ξ−η|1−τ(|ξ|+|η|)p−2+τ.

(2.10)(|ξ|p−2ξ−|η|p−2η)⋅(ξ−η)≥C|ξ−η|2+τ(|ξ|+|η|)p−2−τ.

Lemma 2.2

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

(2.11) |as−bs|≤|a−b|(a+b)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) an+τ∑k=0nbk≤τ∑k=0nγkak+τ∑k=0nck+B,forn≥0.

Suppose that τγ_k < 1, for all k, and set σk=1(1−τγk). Then

(2.13)an+τ∑k=0nbk≤expτ∑k=0nγkσkτ∑k=0nck+B,forn≥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

−Δv=f,inΩ,v=0,on∂Ω,

where Ω⊂Rd, d = 2, 3, is a smooth and bounded domain. Then,

(2.14)‖v‖Hk≤C‖f‖Hk−2,k=2,3.

Let Th={K} be a uniform triangular or tetrahedral partition of Ω into triangles or tetrahedrons. Thus, let h=maxK∈Th{hK} denote the mesh size, where h_K = diam(K) = max{‖x − y‖, x, y ∈ K}, and V_h be the finite dimensional subspace of H01(Ω), which consists of continuous piecewise polynomials of degree r ≥ 1 on Th.

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 {wn}n=0N define

D1wn=wn−wn−1Δt,w¯n=12(wn+wn−1),n=1,2…,N,D2wn=32D1wn−12D1wn−1andw^n=2wn−1−wn−2,n=2,…,N.

The following telescope formula is for n ≥ 2

(2.15)(D2wn,wn)=14Δt‖wn‖2−‖wn−1‖2+‖w^n+1‖2−‖w^n‖2+‖wn−2wn−1+wn−2‖2.

Under the above notations, we propose the following linearized second-order BDF Galerkin finite element scheme associated to (1.1), which is to find Uhn∈Vh such that

• Step 1: For Uh0=Πhu0∈Vh, find Uh1∈Vh such that for all v_h ∈ V_h

(2.16)(D1Uh1,vh)+a(l(U^h1))(∇U¯h1,∇vh)+α(|U^h1|p−2U¯h1,vh)=(f(U^h1),vh),

where U^h1 is given by

(2.17)U^h1−Uh0Δt/2,vh+a(l(Uh0)(∇U^h1,∇vh)+α(|Uh0|p−2U^h1,vh)=(f(Uh0),vh).

• Step 2: For 2 ≤ n ≤ N, find Uhn∈Vh such that for all v_h ∈ V_h

(2.18)(D2Uhn,vh)+a(l(U^hn))(∇Uhn,∇vh)+α(|U^hn|p−2Uhn,vh)=(f(U^hn),vh).

Π_h is an interpolation operator from H01(Ω) 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 Uhn which satisfies

(2.19)‖Uhn‖2+‖U^hn+1‖2+CΔt∑k=1n‖∇Uhk‖2≤C‖u0‖2.

Proof. 1 For the existence, taking vh=Uh1, vh=U^h1 and vh=Uhn in (2.16)–(2.18), respectively, the existence and uniqueness of Uh1, U^h1 and Uhn are from the Lax–Milgram theorem and the hypothesis (H3).

Let vh=U¯hn in (2.16), we have

12Δt(‖Uh1‖2−‖Uh0‖2)+a(l(U^h1))‖∇U¯h1‖2+α(|U^h1|p−2U¯h1,U¯h1)=(f(U^h1),U¯h1)≤‖f(U^h1)‖‖U¯h1‖.

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

12Δt(‖Uh1‖2−‖Uh0‖2)+M‖∇U¯h1‖2≤L‖U^h1‖‖U¯h1‖≤C‖U^h1‖2+M2‖∇U¯h1‖2.

and

(2.20)‖Uh1‖2+CΔt‖∇U¯h1‖2≤CΔt‖U^h1‖2+‖u0‖2.

Now, let vh=U^h1 in (2.17), the same arguments used above give us

‖U^h1‖2+CΔt‖∇U^h1‖2≤C‖u0‖2.

Taking vh=Uhn in (2.18), using the lower bound of a(⋅), (H2) and dropping the third term of the left hand side lead to

(D2Uhn,Uhn)+M‖∇Uhn‖2≤L‖U^hn‖‖Uhn‖.

From the telescope (2.15), we obtain

14Δt(‖Uhn‖2−‖Uhn−1‖2+‖U^hn+1‖2−‖U^hn‖2)+M‖∇Uhn‖2≤CL‖U^hn‖‖∇Uhn‖≤C‖U^hn‖2+M2‖∇Uhn‖2.

That is

(2.21)‖Uhn‖2−‖Uhn−1‖2+‖U^hn+1‖2−‖U^hn‖2+M2Δt‖∇Uhn‖2≤C‖U^hn‖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 Uhn, and

(2.22) max0≤n≤N(‖un−Uhn‖2+‖u^n−U^hn‖2+Δt∑k=0n‖∇(un−Uhn)‖2)≤C(Δt4+h2r+2),

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)D1U1−a(l(U^1))ΔU¯1+α|U^1|p−2U¯1=f(U^1)U1=0on∂Ω,

where U^1 is the solution to

(3.2)U^1−U0Δt/2−a(l(U0))ΔU^1+α|U0|p−2U^1=f(U0)U^1=0on∂Ω.

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

(3.3)D2Un−a(l(U^n))ΔUn+α|U^n|p−2Un=f(U^n)Un=0on∂Ω.

The weak formulations of (3.1)–(3.3) are defined as follows: find Un∈H01(Ω) such that for all v∈H01(Ω)

(3.4)(D1U1,v)+a(l(U^1))(∇U¯1,∇v)+α(|U^1|p−2U¯1,v)=(f(U^1),v),

and for 2 ≤ n ≤ N

(3.5)(D2Un,v)+a(l(U^n))(∇Un,∇v)+α(|U^n|p−2Un,v)=(f(U^n),v),

with U^1∈H01(Ω) such that

(3.6)U^1−U0Δt/2,v+a(l(U0))(∇U^1,∇v)+α(|U0|p−2U^1,v)=(f(U0),v).

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)D1u1−a(l(u^1))Δu¯1+α|u^1|p−2u¯1=f(u^1)+R1

(3.8)D2un−a(l(u^n))Δun+α|u^n|p−2un=f(u^n)+Rn,n=2,…,N,

where u^1 satisfies

(3.9)u^1−u0Δt/2−a(l(u0))Δu^1+α|u0|p−2u^1=f(u0)+R^0.

R^0, R¹ and Rⁿ are, respectively, the truncation errors given by

R^0=u^1−u0Δt/2−ut1/2−Δu^1(a(l(u0))−a(l(u1/2)))+α|u0|p−2(u1/2−u0)+α(|u0|p−2u0−|u1/2|p−2u1/2)+(f(u0)−f(u1/2)),R1=(D11u1−ut1/2)−Δu¯1(a(l(u^1))−a(l(u1/2)))−a(l(u1/2))Δ(u¯1−u1/2)+α|u^1|p−2(u¯1−u^1)+α(|u^1|p−2u^1−|u1/2|p−2u1/2)+(f(u^1)−f(u1/2))Rn=(D2un−utn)−Δun(a(l(u^n))−a(l(un)))+α|u^n|p−2(un−u^n)+α(|u^n|p−2u^n−|un|p−2un)+(f(u^n)−f(un)).

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

(3.10)∑n=1NΔt‖Rn‖121/2+Δt‖R^0‖1≤CΔt2.

Let us denote

e¯1=u¯1−U¯1;e^n=u^n−U^n;en=un−Unfor1≤n≤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)‖e^1‖+Δt1/2‖∇e^1‖≤CΔt2.

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

2e^1−Δta(l(u0))Δe^1+αΔt|u0|p−2e^1=ΔtR^0.

Testing the above equation by e^1 yield

2‖e^1‖2+Δta(l(u0))‖∇e^1‖2+αΔt∫Ω|u0|p−2|e^1|2dx=Δt(R^0,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‖e^1‖2+mΔt‖∇e^1‖2+αΔt∫Ω|u0|p−2|e^1|2dx≤12Δt2‖R^0‖2+12‖e^1‖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) ‖e1‖2+Δt‖∇e1‖2≤CΔt5.

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

e1−Δta(l(u0))Δe¯1+αΔt|u^1|p−2e¯1=−αΔtu¯1(|u^1|p−2−|U^1|p−2)+ΔtΔu¯1(a(l(u^1))−a(l(U^1)))+Δt(f(u^1)−f(U^1))+ΔtR1.

Testing the above equation by e¹ and using the fact that e¯1=12e1, we have

(3.13)‖e1‖2+12mΔt‖∇e1‖2+12αΔt∫Ω|u^1|p−2|e1|2dx≤αΔt(u¯1(|u^1|p−2−|U^1|p−2),e1)+Δt(a(l(u^1))−a(l(U^1)))(∇u¯1,∇e1)+Δt(f(u^1)−f(U^1),e1)+Δt(R1,e1)=∑i=14Ii.

We have

I1=αΔt(u¯1(|u^1|p−2−|U^1|p−2),e1)≤C1(‖u¯1‖L∞,‖u^1‖L∞,‖U^1‖L∞,p)Δt‖e^1‖‖e1‖≤C1Δt2‖e^1‖2+14‖e1‖2≤CΔt6+14‖e1‖2(using (3.11)).

I2=Δt(a(l(u^1))−a(l(U^1)))(∇u¯1,∇e1)≤ΔtA‖e^1‖‖∇u¯1‖‖∇e1‖(using (H4))≤C2(A,‖∇u¯1‖)Δt‖e^1‖2+m4Δt‖∇e1‖2≤CΔt5+m4Δt‖∇e1‖2.

I3=Δt(f(u^1)−f(U^1),e1)≤γΔt‖e^1‖‖e1‖(using (H2))≤C3Δt2‖e^1‖2+14‖e1‖2≤CΔt6+14‖e1‖2.

I4=Δt(R1,e1)≤C4‖R1‖2+14‖e1‖2≤CΔt5+14‖e1‖2(using (3.11)).

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) max1≤n≤N‖en‖2+‖e^n‖2+Δt∑k=1n‖∇ek‖2≤CΔt4,

(3.15)max0≤n≤N‖Un‖∞≤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 U^n and the induction assumption, ‖U^n‖∞≤C.

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

(3.16)D2en−a(l(U^n))Δen+α|U^n|p−2en=(a(l(u^n))−a(l(U^n)))Δun−αun(|u^n|p−2−|U^n|p−2)+(f(u^n)−f(U^n))+Rn.

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

‖en‖2−‖en−1‖2+‖e^n+1‖2−‖e^n‖2+‖en−2en−1+en−2‖2+4Δta(l(U^n))‖∇en‖2+4Δtα∫Ω|U^n|p−2|en|2dx=−4Δt(a(l(u^n))−a(l(U^n)))(∇un,∇en)−4Δtα(un(|u^n|p−2−|U^n|p−2),en)+4Δt(f(u^n)−f(U^n),en)+4Δt(Rn,en).

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

(3.17)‖en‖2−‖en−1‖2+‖e^n+1‖2−‖e^n‖2+4mΔt‖∇en‖2≤4Δt(a(l(u^n))−a(l(U^n)))(∇un,∇en)+4Δtα(un(|u^n|p−2−|U^n|p−2),en)+4Δt(f(u^n)−f(U^n),en)+4Δt(Rn,en)=∑k=14Jk.

We have

J1=4Δt(a(l(u^n))−a(l(U^n)))(∇un,∇en)≤4ΔtA‖e^n‖‖∇un‖‖∇en‖(using (H4))≤CΔt‖e^n‖2+mΔt‖∇en‖2

J2=4Δtα(un(|u^n|p−2−|U^n|p−2),en)≤C(‖un‖L∞,‖u^n‖L∞,‖U^n‖L∞,p)Δt‖e^n‖‖en‖≤CΔt‖e^n‖2+‖en‖2

J3=4Δt(f(u^n)−f(U^n),en)≤4γΔt‖e^n‖‖en‖(using (H2))≤CΔt‖e^n‖2+‖en‖2.

J4=4Δt(Rn,en)≤CΔt‖Rn‖2+‖en‖2.

Therefore,

‖en‖2−‖en−1‖2+‖e^n+1‖2−‖e^n‖2+Δt‖∇en‖2≤CΔt‖en‖2+‖e^n‖2+‖Rn‖2.

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

‖en‖2+‖e^n‖2+Δt∑k=1n‖∇ek‖2≤CΔt4.

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

‖Un‖≤‖un‖+‖en‖≤C,‖DkUn‖≤‖Dkun‖+‖Dken‖≤C,withk=1,2.

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

‖Un‖2≤C‖DkUn‖+C‖∇Un‖+C‖|U^n|p−2Un‖+L‖U^n‖≤C+C‖U^n‖∞p−2‖Un‖≤C.

Using (2.3), we have

‖Un‖∞≤C‖Un‖1/2‖Un‖21/2≤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 Rh:H01(Ω)→Vh be a Ritz projection operator defined by

(4.1)∫Ω∇(u−Rhu)⋅∇wdx=0,∀w∈H01(Ω).

Then we have the following lemma.

Lemma 4.1

(cf. Ref. [19]). If u∈Hr+1(Ω)∩H01(Ω), then

(4.2) ‖u−Πhu‖+h‖∇(u−Πhu)‖≤Chr+1‖u‖Hr+1,

(4.3) ‖u−Rhu‖+h‖∇(u−Rhu)‖≤Chr+1‖u‖Hr+1

(4.4) ‖vh‖∞≤Ch−d/2‖vh‖,∀vh∈Vh.

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

Let us denote

Eh0=u0−Uh0e¯h1=RhU¯1−U¯h1;e^hn=RhU^n−U^hn;ehn=RhUn−UhnE¯1=RhU¯1−U¯1;E^n=RhU^n−U^n;En=RhUn−Unfor1≤n≤N.

From lemma 4.1, we have

(4.5)‖Eh0‖+h‖∇Eh0‖≤Chr+1

(4.6)‖E¯1‖+‖E^n‖+‖En‖+h(‖∇E¯1‖+‖∇E^n‖+‖∇En‖)≤Chr+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) ‖e^h1‖2+Δt‖∇e^h1‖2≤Ch2r+2

(4.8)‖U^h1‖≤C.

Proof. From equations (2.17) and (3.6), e^h1 satisfies the following equation:

2(e^h1,vh)+Δta(l((Uh0))(∇e^h1,∇vh)+αΔt(|Uh0|p−2e^h1,vh)=2(Eh0,vh)−2(E^1,vh)−Δt(a(l((U0))−a(l((Uh0)))(∇U^1,∇vh)−αΔt(|Uh0|p−2E^1,vh)−αΔt((|U0|p−2−|Uh0|p−2)U^1,vh)+Δt(f(U0)−f(Uh0),vh).

Setting vh=e^h1 in the above equations leads to

(4.9)2‖e^h1‖2+mΔt‖∇e^h1‖2+αΔt∫Ω|Uh0|p−2|e^h1|2dx≤2(Eh0+E^1,e^h1)+Δt(a(l((U0))−a(l((Uh0)))(∇U^1,∇e^h1)+αΔt(|Uh0|p−2E^1,e^h1)+αΔt((|U0|p−2−|Uh0|p−2)U^1,e^h1)+Δt(f(U0)−f(Uh0),e^h1)=∑i=15Li.

From (4.5) and (4.6), we have

L1=2(Eh0+E^1,e^h1)≤Ch2r+2+12‖e^h1‖2,L2=Δt(a(l((U0))−a(l((Uh0)))(∇U^1,∇e^h1)≤ΔtA‖Eh0‖‖∇U^1‖‖∇e^h1‖≤CΔth2r+2+m2Δt‖∇e^h1‖2,L3=αΔt(|Uh0|p−2E^1,e^h1)≤αΔtC(‖Uh0‖∞)‖E^1‖‖e^h1‖≤CΔt2h2r+2+14‖e^h1‖2,L4=αΔt((|U0|p−2−|Uh0|p−2)U^1,e^h1)≤C(‖U0‖∞,‖Uh0‖∞,‖U^1‖∞)Δt‖Eh0‖‖e^h1‖≤CΔt2h2r+2+14‖e^h1‖2,L5=Δt(f(U0)−f(Uh0),e^h1)≤γΔt‖Eh0‖‖e^h1‖≤CΔt2h2r+2+14‖e^h1‖2.

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

From the inverse inequality, ‖e^h1‖∞≤Ch−d/2‖e^h1‖≤Chr/2 and

‖U^h1‖∞≤‖RhU1‖∞+‖e^h1‖∞≤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) max1≤n≤N‖ehn‖2+‖e^hn‖2+Δt∑k=1n‖∇ehk‖2≤Ch2r+2,

(4.11) max0≤n≤N‖Uhn‖≤C.

Proof. The proof of this result will be done by mathematical induction. Since Uh0=ΠhU0, (4.11) holds for n = 0. To compute the error estimate (4.10) for n = 1, subtract (3.4) from (2.16) and take vh=e¯h1≡12(eh1+Eh0),

(4.12)‖eh1‖2+2mΔt‖∇e¯h1‖2+2αΔt∫Ω|U^h1|p−2|e¯h1|2dx≤(E1+Eh0,e¯h1)+2Δt(a(l((U^1))−a(l((U^h1)))(∇U¯1,∇e¯h1)+2Δtα(|U^h1|p−2E¯1,e¯h1)+2Δtα((|U^1|p−2−|U^h1|p−2)U¯1,e¯h1)+2Δt(f(U^1)−f(U^h1),e¯h1)=∑i=15Ti.

From (4.5) and (4.6), we have

T1=(E1+Eh0,e¯h1)≤Ch2r+2+18‖eh1‖2,T2=2Δt(a(l((U^1))−a(l((U^h1)))(∇U¯1,∇e¯h1)≤ΔtA‖U^1−U^h1‖‖∇U¯1‖‖∇e¯h1‖≤CΔth2r+2+CΔt‖e^h1‖2+m2Δt‖∇e¯h1‖2,T3=2αΔt(|U^h1|p−2E¯1,e¯h1)≤αΔtC(‖U^h1‖∞)‖E¯1‖‖e¯h1‖≤CΔt2h2r+2+18‖eh1‖2,T4=2αΔt((|U^1|p−2−|U^h1|p−2)U¯1,e¯h1)≤C(‖U^1‖∞,‖U^h1‖∞,‖U¯1‖∞)Δt‖U^1−U^h1‖‖e¯h1‖≤CΔt2h2r+2+CΔt2‖e^h1‖2+18‖eh1‖2,T5=2Δt(f(U^1)−f(U^h1),e¯h1)≤γΔt‖U^1−U^h1‖‖e¯h1‖≤CΔt2h2r+2+CΔt2‖e^h1‖2+18‖eh1‖2.

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

(4.13)‖eh1‖2+Δt‖∇e¯h1‖2≤Ch2r+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 U^hn and the induction assumption, ‖U^hn‖∞≤C.

Subtracting (2.18) from (3.5), we obtain

(D2ehn,vh)+a(l((e^hn))(∇ehn,∇vh)+α(|e^hn|p−2ehn,vh)=−(D2En,vh)−(a(l((U^n))−a(l((U^hn)))(∇Un,∇vh)−α(|U^hn|p−2En,vh)−α((|U^n|p−2−|U^hn|p−2)Un,vh)+(f(U^n)−f(U^hn),vh).

If one takes vh=4Δtehn and uses the telescope formula, one obtains

‖ehn‖2−‖ehn−1‖2+‖e^hn+1‖2−‖e^hn‖2+‖ehn−2ehn−1+ehn−2‖2+4Δta(l(U^hn))‖∇ehn‖2+4Δtα∫Ω|U^hn|p−2|ehn|2dx=−4Δt(D2En,ehn)−4Δt(a(l(U^n))−a(l(U^hn)))(∇Un,∇ehn)−4Δtα(|U^hn|p−2En,ehn)−4Δtα((|U^n|p−2−|U^hn|p−2)Un,ehn)+4Δt(f(U^n)−f(U^hn),ehn).

That is

(4.14)‖ehn‖2−‖ehn−1‖2+‖e^hn+1‖2−‖e^hn‖2+4mΔt‖∇ehn‖2≤4Δt(D2En,ehn)+4Δt(a(l(U^n))−a(l(U^hn)))(∇Un,∇ehn)+4Δtα(|U^hn|p−2En,ehn)+4Δtα((|U^n|p−2−|U^hn|p−2)Un,ehn)+4Δt(f(U^n)−f(U^hn),ehn)=∑i=15Ki.

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

K1=4Δt(D2En,ehn)≤CΔt(‖D2En‖2+‖ehn‖2)≤CΔth2r+2+CΔt‖ehn‖2,K2=4Δt(a(l((U^n))−a(l((U^hn)))(∇Un,∇ehn)≤ΔtA‖E^n+e^hn‖‖∇Un‖‖∇ehn‖≤CΔth2r+2+CΔt‖e^hn‖2+mΔt‖∇ehn‖2,K3=4αΔt(|U^hn|p−2En,ehn)≤4αΔtC(‖U^hn‖∞)‖En‖‖ehn‖≤CΔth2r+2+CΔt‖ehn‖2,K4=4αΔt((|U^n|p−2−|U^hn|p−2)Un,ehn)≤C(‖U^n‖∞,‖U^hn‖∞,‖Un‖∞)Δt‖E^n+e^hn‖‖ehn‖≤CΔth2r+2+CΔt(‖e^hn‖2+‖ehn‖2),K5=4Δt(f(U^n)−f(U^hn),ehn)≤4γΔt‖E^n+e^hn‖‖ehn‖≤CΔth2r+2+CΔt(‖e^hn‖2+‖ehn‖2).

Taking these bounds into (4.14), we obtain

‖ehn‖2−‖ehn−1‖2+‖e^hn+1‖2−‖e^hn‖2+4mΔt‖∇ehn‖2≤CΔth2r+2+CΔt(‖e^hn‖2+‖ehn‖2).

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

‖ehn‖2+‖e^hn‖2+Δt∑k=1n‖∇ehk‖2≤exp(CT)h2r+2.

From the inverse inequality, ‖ehn‖∞≤Ch−d/2‖ehn‖≤Chr/2 and

‖Uhn‖∞≤‖RhUn‖∞+‖ehn‖∞≤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)ut−a(l(u))Δu+α|u|p−2u=f(u)+ginΩ×(0,T]u(x,t)=0on∂Ω×(0,T]u(x,0)=u0(x)inΩ,

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

f(s)=s(10−s),a(s)=3+cos(s),l(u)=∫Ωudx.

g is chosen correspondingly to the exact solution

u(x,y,t)=2(1+t2exp(−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