A Galerkin Finite Element Method for a Nonlocal Parabolic System with Nonlinear Boundary Conditions Arising from the Thermal Explosion Theory

Abstract

In this paper, we discuss a class of nonlocal parabolic systems with nonlinear boundary conditions arising from the thermal explosion theory. First, we prove the local existence and uniqueness of the classical solution using the Leray–Schauder fixed-point theorem. Then, we analyze three Galerkin approximations of the system and derive the optimal-order error estimates: $\mathcal{O}(h^{r+1})$ in $L^2$ norm for continuous-time Galerkin approximation, $\mathcal{O}(h^{r+1}+{(\Delta t)}^2)$ in the $L^2$ norm for Crank–Nicolson Galerkin approximation, and $\mathcal{O}(h^{r+1}+{(\Delta t)}^2)$ in both $L^2$ and $H^1$ norms for extrapolated Crank–Nicolson Galerkin approximation.

1. Introduction

In this paper, we discuss the following semi-linear parabolic equations:

(1)$\begin{array}{c}\left\{\begin{array}{cc}{\displaystyle \frac{\partial u}{\partial t}}=\Delta u+k_1(x,t){\int }_{\Omega }{f}_1(u(x,t),v(x,t))\mathrm{d}x,\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ {\displaystyle \frac{\partial v}{\partial t}}=\Delta v+k_2(x,t){\int }_{\Omega }{f}_2(u(x,t),v(x,t))\mathrm{d}x,\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ {\displaystyle \frac{\partial u}{\partial \mathit{n}}}+g_1(u)u={\displaystyle \frac{\partial v}{\partial \mathit{n}}}+g_2(v)v=0,\hfill & (x,t)\in \partial \Omega \times (0,T],\hfill \\ u(x,0)=u_0(x),v(x,0)=v_0(x),\hfill & x\in \Omega ,\hfill \end{array}\right.\hfill \end{array}$

where $\Omega$ is a smooth, convex, bounded domain in ${\mathbb{R}}^d(d\ge 1)$ with a smooth boundary $\partial \Omega$, ${\displaystyle \frac{\partial u}{\partial \mathit{n}}}$ and ${\displaystyle \frac{\partial v}{\partial \mathit{n}}}$ are the outer normal directional derivatives of u and v on $\partial \Omega$, respectively, and $k_1$, $k_2$, $f_1$, $f_2$, $g_1$, $g_2$, $u_0$, and $v_0$ are known functions satisfying some assumptions.

Motivated by [1,2,3,4], we consider the nonlocal parabolic system with nonlinear boundary conditions in the thermal explosion theory. As noted in [4], in certain thermal explosion problems that involve prolonged induction times (such as the safe storage of energetic materials or nuclear waste), the conventional Dirichlet boundary condition $u=v=0$ is no longer valid. This is because, during this period, the temperature of the reactive material is significantly higher than that of the surrounding environment. Consequently, it is necessary to apply heat loss boundary conditions, as illustrated in Equation (1), to accurately describe the temperature distribution at the boundary. Our focus here is on the numerical solution of system, where we employ the Galerkin method to approach the problem.

The Galerkin method is a widely used numerical technique for solving partial differential equations. It involves using the weak form of the original equation, followed by subdividing the region into smaller elements using piecewise polynomials in the finite element approximation space. Polynomials are then used to approximate the unknown functions on each element, and, ultimately, solvable linear equations are derived. For example, in [5], some results on elliptic equations are presented; linear second-order hyperbolic equations with Dirichlet boundary conditions are discussed in [6], and in [7], nonlinear hyperbolic equations with non-homogeneous boundary conditions are studied, with their superconvergence being further analyzed in [8].

The parabolic equation also has a lot of results. In [9], the authors analyze the following semilinear parabolic equations for homogeneous boundary conditions:

${\displaystyle \frac{\partial u_j(x,t)}{\partial t}}-{\alpha }_j\nabla u_j(x,t)+R_j(u(x,t),x,t)=0,j=1,\cdots ,M.$

In [10], the nonlinear parabolic equation is extended to the following form:

$c(x,u){\displaystyle \frac{\partial u}{\partial t}}-\nabla ·a(x,u)\nabla u+\sum _{i=1}^p{b}_i(x,u){\displaystyle \frac{\partial u}{\partial x_i}}+f(x,u)=0.$

In [11,12], using a nonlinear elliptic projection, the following linear and nonlinear parabolic equations are treated:

${\displaystyle \frac{\partial u}{\partial t}}=\nabla ·(a\nabla u),$

$c(x,u){\displaystyle \frac{\partial u}{\partial t}}-\nabla ·(a(x,u)\nabla u)=f(x,t,u)$

with nonlinear boundary conditions $a(x){\displaystyle \frac{\partial u}{\partial \mathit{n}}}=g\left(x,t,u\right)$, $a\left(x,u\right){\displaystyle \frac{\partial u}{\partial \mathit{n}}}=g\left(x,t,u\right)$, respectively. In [13], the author modified the elliptic projection proposed in [11] and successfully applied it to parabolic integro-differential equations with nonlinear boundary conditions

${\displaystyle \frac{\partial u}{\partial t}}=\nabla ·\left\{a(u)\nabla u+{\int }_0^{t}b(t,u(x,s))\nabla u(x,s)ds\right\}+f(u).$

In recent years, the Galerkin finite element method has also been applied to various types of equations, such as fractional partial differential equations (see [14,15]), stochastic differential equations (see [16,17]), etc. At the same time, it has also been continuously developed as the discontinuous Galerkin method [18,19], $H^1$-Galerkin mixed finite element method [20,21], and spectral Galerkin method [22,23]. In addition, there are some interesting results related to our work that can be found in [24,25,26,27,28,29].

In this paper, we apply the Galerkin method to nonlocal parabolic systems with nonlinear boundary conditions for the first time, thus extending the equation form in [12] to the case with nonlinear nonlocal heat sources. Three Galerkin approximations were successfully proposed, and the existence, uniqueness and optimal-order error estimates for each of these approximations were obtained.

Before discussing the approximate solution, it is necessary to confirm the existence and uniqueness of the classical solution. However, there is currently no suitable reference, so we provide a detailed demonstration in a separate section. Before discussing the classical solution of System (1), we provide the following hypothesis:

(A1) $f_i,g_i^{\prime }(i=1,2)$ satisfy the local Lipschitz condition;

(A2) $k_i(i=1,2)$ satisfy the Holder condition $\left|k_i(x,t)-k_i(y,s)\right|\le C\left({|x-y|}^{\alpha }+{|t-s|}^{{\displaystyle \frac{\alpha }{2}}}\right)$, where $(0<\alpha <1)$;

(A3) $u_0,v_0\in C^{2+\alpha }(\overline{\Omega })(0<\alpha <1)$ are non-negative.

This article is arranged as follows: In Section 1, we introduce the necessary notation and useful lemmas. In Section 2, we prove the local existence and uniqueness of the classical solution using the Leray–Schauder fixed-point theorem. Section 3, Section 4 and Section 5 present the continuous-time Galerkin approximation, Crank–Nicolson Galerkin approximation, and extrapolated Crank–Nicolson Galerkin approximation, respectively. The uniqueness and optimal error estimation of the three numerical schemes are also derived.

In the following discussion, we set

$\begin{array}{cc}(f,g)={\int }_{\Omega }fg\mathrm{d}x,\hfill & {\parallel f\parallel }^2=(f,f),\hfill \\ ⟨f,g⟩={\int }_{\partial \Omega }fg\mathrm{d}x,\hfill & {|f|}^2=⟨f,f⟩.\hfill \end{array}$

Throughout this paper, we use the standard notation $W^{k,p}(\Omega )$ for a Sobolev space on $\Omega$. If $p=2$, we usually write $H^k(\Omega )=W^{k,2}(\Omega )$, and we donate the norms on $\Omega$ and $\partial \Omega$ via ${∥·∥}_{H^k\left(\Omega \right)}$ and ${\left|·\right|}_{H^k\left(\partial \Omega \right)}$, respectively. Let X be a Banach space and $\varphi (t):[0,T]\mapsto X$. We define the following norms:

$\begin{array}{cc}\hfill & {\parallel \varphi \parallel }_{L^p(a,b,X)}={\left[{\int }_a^b{\parallel \varphi (s)\parallel }_X^p\mathrm{d}s\right]}^{{\displaystyle \frac{1}{p}}},1\le p<\infty ,\hfill \\ \hfill & {\parallel \varphi \parallel }_{L^{\infty }(a,b,X)}=\underset{a\le s\le b}{\mathrm{ess}\sup }{\parallel \varphi (s)\parallel }_X.\hfill \end{array}$

When $a=0,b=T$, these are often abbreviated as ${∥\varphi ∥}_{L^p\left(X\right)}$ and ${∥\varphi ∥}_{L^{\infty }\left(X\right)}$, respectively.

In the following sections, we have several commonly used results.

Lemma 1 

([30], p. 258, Trace Theorem).

$\begin{array}{c}{|\varphi |}^2\le {\epsilon \parallel \nabla \varphi \parallel }^2+C(\epsilon ){\parallel \varphi \parallel }^2(\forall \epsilon >0),\\ {|\varphi |}^2\le C\parallel \nabla \varphi \parallel \parallel \varphi \parallel .\end{array}$

Lemma 2. 

${∥{\int }_{\Omega }\varphi \left(x,t\right)\mathrm{d}x∥}^2\le C{∥\varphi \left(x,t\right)∥}^2$

Proof of Lemma 2. 

The conclusion can be obtained by using the definition of the $L^2$ norm and Cauchy–Schwarz inequality. □

2. The Existence and Uniqueness of the Classical Solution

In this section, we establish the existence and uniqueness of the classical solution to System (1). Although this result is an application of the fixed-point theorem, a comprehensive proof for this specific problem has not been previously provided. Therefore, we present a detailed discussion here, which is a modification of the proof of Theorem 1.1 in reference [31]. To begin, we introduce several symbols that will be used in this section (see [32]). For $l\in \mathbb{N}$ and $\alpha \in (0,1]$, we define the following norms and seminorms:

$\begin{array}{cc}\hfill {\left[f\right]}_{l+\alpha ,{\Omega }_T}& :=\sum _{|\beta |+2j=l}\underset{(x,t)\ne (y,s)\in {\Omega }_T}{sup}{\displaystyle \frac{\left|{\partial }_x^{\beta }{\partial }_t^j(f(x,t)-f(y,s))\right|}{{|x-y|}^{\alpha }+{|t-s|}^{\frac{\alpha }{2}}}},\hfill \\ \hfill {⟨f⟩}_{l+\alpha ,{\Omega }_T}& :=\sum _{|\beta |+2j=l-1}\underset{(x,t)\ne (x,s)\in {\Omega }_T}{sup}{\displaystyle \frac{\left|{\partial }_x^{\beta }{\partial }_t^j(f(x,t)-f(x,s))\right|}{{|t-s|}^{\frac{1+\alpha }{2}}}},\hfill \end{array}$

$\begin{array}{cc}\hfill {|f|}_{l+\alpha ,{\Omega }_T}& :={\left[f\right]}_{l+\alpha ,{\Omega }_T}+{⟨f⟩}_{l+\alpha ,{\Omega }_T},\hfill \\ \hfill {\parallel f\parallel }_{C^l\left({\Omega }_T\right)}& :=\sum _{|\beta |+2j\le l}\underset{(x,t)\in {\Omega }_T}{sup}\left|{\partial }_x^{\beta }{\partial }_t^{j}f(x,t)\right|,\hfill \\ \hfill {\parallel f\parallel }_{C^{l+\alpha }\left({\Omega }_T\right)}& :={\parallel f\parallel }_{C^l\left({\Omega }_T\right)}+{|f|}_{l+\alpha ,{\Omega }_T},\hfill \end{array}$

where ${\Omega }_T:=\Omega \times (0,T]$. Now, we define the following functional spaces:

$\begin{array}{cc}\hfill C^l\left({\Omega }_T\right)& :=\left\{f:{\partial }_x^{\beta }{\partial }_t^{j}f\mathrm{is}\mathrm{continuous}\mathrm{in}{\Omega }_T\mathrm{for}|\beta |+2j=l\right\},\hfill \\ \hfill C^{l+\alpha }\left({\Omega }_T\right)& :=\left\{f\in C^l\left({\Omega }_T\right):{\parallel f\parallel }_{C^{l+\alpha }\left({\Omega }_T\right)}<\infty \right\}.\hfill \end{array}$

One can verify that $C^l\left({\overline{\Omega }}_T\right)$ and $C^{l+\alpha }\left({\overline{\Omega }}_T\right)$ are Banach spaces. The boundary smoothness condition is necessary to guarantee the inclusion $C^a\left({\overline{\Omega }}_T\right)\subset C^b\left({\overline{\Omega }}_T\right)$ for $b>a\ge 0$, since this is not generally true for an arbitrary domain $\Omega$ (see [33], p. 53). Moreover, $C^a\left({\overline{\Omega }}_T\right)$ is compactly embedded in $C^b\left({\overline{\Omega }}_T\right)$ for any $a\ge 1$ and $a>b\ge 0$ (see [33], Lemma 6.36).

Now, let us recall a useful result for the linear model. We consider the linear second-order parabolic equation of the non-divergence form as follows:

(2)$\left\{\begin{array}{c}{\mathcal{L}}_{1}u:={\displaystyle \frac{\partial u}{\partial t}}-\sum _{i,j=1}^n{a}_{ij}(x,t)D_{ij}u+\sum _{i=1}^n{b}_i(x,t)D_{i}u+c(x,t)u=F_1(x,t),\mathrm{in}{\Omega }_T,\hfill \\ {\mathcal{L}}_{2}v:={\displaystyle \frac{\partial v}{\partial t}}-\sum _{i,j=1}^n{\widehat{a}}_{ij}(x,t)D_{ij}v+\sum _{i=1}^n{\widehat{b}}_i(x,t)D_{i}v+\widehat{c}(x,t)v=F_2(x,t),\mathrm{in}{\Omega }_T.\hfill \end{array}\right.$

Assume that there exists $\Lambda \ge \lambda >0$, $\widehat{\Lambda }\ge \widehat{\lambda }>0$, such that

(3)${\lambda |\xi |}^2\le \sum _{i,j=1}^n{a}_{ij}(x,t){\xi }_i{\xi }_j\le \Lambda {|\xi |}^2,(x,t)\in {\Omega }_T,\xi \in {\mathbb{R}}^n,$

(4)$\widehat{\lambda }{|\xi |}^2\le \sum _{i,j=1}^n{\widehat{a}}_{ij}(x,t){\xi }_i{\xi }_j\le \widehat{\Lambda }{|\xi |}^2,(x,t)\in {\Omega }_T,\xi \in {\mathbb{R}}^n,$

where $a_{ij},{\widehat{a}}_{ij},b_i,{\widehat{b}}_i,c,\widehat{c}\in C^{\alpha }\left({\overline{\Omega }}_T\right)(0<\alpha <1)$ and

(5)${\displaystyle \frac{1}{\lambda }}\left\{\sum _{i,j=1}^n{∥a_{ij}∥}_{C^{\alpha }\left({\overline{\Omega }}_T\right)}+\sum _{i=1}^n{∥b_i∥}_{C^{\alpha }\left({\overline{\Omega }}_T\right)}+{\parallel c\parallel }_{C^{\alpha }\left({\overline{\Omega }}_T\right)}\right\}\le {\Lambda }_{\alpha },$

(6)${\displaystyle \frac{1}{\widehat{\lambda }}}\left\{\sum _{i,j=1}^n{∥{\widehat{a}}_{ij}∥}_{C^{\alpha }\left({\overline{\Omega }}_T\right)}+\sum _{i=1}^n{∥{\widehat{b}}_i∥}_{C^{\alpha }\left({\overline{\Omega }}_T\right)}+{\parallel \widehat{c}\parallel }_{C^{\alpha }\left({\overline{\Omega }}_T\right)}\right\}\le {\widehat{\Lambda }}_{\alpha }.$

Theorem 1 

([32], p. 79, Theorem 4.31). Let Assumptions (3)–(6) be in force, and $\partial \Omega \in$$C^{2+\alpha })$$(0<\alpha <1$. Let $F_1,F_2\in C^{\alpha }\left({\overline{\Omega }}_T\right),G_1,G_2\in C^{1+\alpha }\left({\overline{\Omega }}_T\right)$ and $u_0,v_0\in C^{2+\alpha }(\overline{\Omega })$ satisfy the first-order compatibility condition:

(7)${\displaystyle \frac{\partial u_0}{\partial n}}=G_1(x,0),{\displaystyle \frac{\partial v_0}{\partial n}}=G_2(x,0)\mathrm{on}\partial \Omega .$

Then, there exists a unique solution $u,v\in C^{2+\alpha }\left({\overline{\Omega }}_T\right)$ to Problem (28) with the Neumann boundary condition ${\displaystyle \frac{\partial u}{\partial n}}=G_1$, ${\displaystyle \frac{\partial v}{\partial n}}=G_2$ on $\partial \Omega \times (0,T]$. Moreover, there exists a constant C independent of $G_1,G_2$ and $u_0,v_0$, such that

(8)$\begin{array}{cc}\hfill {\parallel (u,v)\parallel }_{C^{2+\alpha }\left({\overline{\Omega }}_T\right)\times C^{2+\alpha }\left({\overline{\Omega }}_T\right)}& \le C\left({\displaystyle \frac{1}{\lambda }}{∥\left(F_1,F_2\right)∥}_{C^{\alpha }\left({\overline{\Omega }}_T\right)\times C^{\alpha }\left({\overline{\Omega }}_T\right)}+{∥\left(G_1,G_2\right)∥}_{C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)}\right)\hfill \\ & +C{∥\left(u_0,v_0\right)∥}_{C^{2+\alpha }(\overline{\Omega })\times C^{2+\alpha }(\overline{\Omega })},\hfill \end{array}$

where C is dependent only on $n,\alpha ,\Lambda /\lambda ,\widehat{\Lambda }/\widehat{\lambda },{\Lambda }_{\alpha },{\widehat{\Lambda }}_{\alpha }$ and Ω.

This estimate, together with the Leray–Schauder fixed-point argument, is the main tool used to prove the following theorem:

Theorem 2. 

Assume that Ω is an open convex bounded domain in ${\mathbb{R}}^d$, where $d\ge 1$ with smooth boundary $\partial \Omega \in C^{\infty }$ and non-negative functions $u_0$, $v_0$ are in $C^{2+\alpha }(\overline{\Omega })$ such that

(9) ${\displaystyle \frac{\partial u_0}{\partial \mathit{n}}}=-g_1(u_0)u_0,{\displaystyle \frac{\partial v_0}{\partial \mathit{n}}}=-g_2(v_0)v_0\mathrm{on}\partial \Omega ,$

where $\alpha \in (0,1)$. Then, there exists $T>0$ such that Problem (1) under the boundary condition

${\displaystyle \frac{\partial u}{\partial \mathit{n}}}=-g_1(u)u,{\displaystyle \frac{\partial v}{\partial \mathit{n}}}=-g_2(v)v\mathrm{on}\partial \Omega \times (0,T],$

admits a unique solution $(u,v)$ in $C^{2+\alpha }\left({\overline{\Omega }}_T\right)\times C^{2+\alpha }\left({\overline{\Omega }}_T\right)$.

Proof. 

From now until the end of this proof, we use C to represent a general constant, which is different in each formula. First, we will prove the local existence of a classical solution using a fixed-point argument. Let $u\in C^{1+\alpha }\left({\overline{\Omega }}_T\right)$ be such that $u(x,0)=u_0(x)$ in $\Omega$ and $v\in C^{1+\alpha }\left({\overline{\Omega }}_T\right)$ be such that $v(x,0)=v_0(x)$ in $\Omega$. Then, the functions $u_0$ and $G_1(x,t)=-g_1(u(x,t))u(x,t)$, $v_0$ and $G_2(x,t)=-g_2(u(x,t))v(x,t)$ satisfy Condition (7), and we can verify $G_1,G_2\in C^{1+\alpha }\left({\overline{\Omega }}_T\right)$.

Assume $T<1$, and consider the set of functions given by

$B_T(M):=\left\{(u,v)\in C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)\mathrm{such}\mathrm{that}max\left\{{\parallel u\parallel }_{C^{1+\alpha }\left({\overline{\Omega }}_T\right)},{\parallel v\parallel }_{C^{1+\alpha }\left({\overline{\Omega }}_T\right)}\right\}\le M\right\}.$

Now, we define the map

$A:B_T(M)⟶(C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)),$

where $A(u,v):=(A_1(u,v),A_2(u,v))=(U,V)$ is a solution of

(10)$\begin{array}{c}\left\{\begin{array}{cc}{U}_t=\Delta U+k_1(x,t){\int }_{\Omega }{f}_1(u(x,t),v(x,t))\mathrm{d}x,\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ V_t=\Delta V+k_2(x,t){\int }_{\Omega }{f}_2(u(x,t),v(x,t))\mathrm{d}x,\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ {\displaystyle \frac{\partial U}{\partial \mathit{n}}}+g_1(u)u=0,{\displaystyle \frac{\partial V}{\partial \mathit{n}}}+g_2(v)v=0,\hfill & (x,t)\in \partial \Omega \times (0,T],\hfill \\ U(x,0)=u_0(x),V(x,0)=v_0(x),\hfill & x\in \Omega .\hfill \end{array}\right.\hfill \end{array}$

We first prove that A sends bounded sets into relative compact sets of $C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)$. Indeed, Inequality (8) implies that there exists $M^{\prime }>0$ independent of T such that

${\parallel A(u,v)\parallel }_{C^{2+\alpha }\left({\overline{\Omega }}_T\right)\times C^{2+\alpha }\left({\overline{\Omega }}_T\right)}\le M^{\prime }$

for all $(u,v)$ in $B_T(M)$. As bounded sets in $C^{2+\alpha }\left({\overline{\Omega }}_T\right)\times C^{2+\alpha }\left({\overline{\Omega }}_T\right)$ are relatively compact in $C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)$. We claim that A is continuous. In fact, let $u_n\to u$ in $C^{1+\alpha }\left({\overline{\Omega }}_T\right)$ and $v_n\to v$ in $C^{1+\alpha }\left({\overline{\Omega }}_T\right)$. Thus, we need to prove $(U_n,V_n):=A(u_n,v_n)\to (U,V):=A(u,v)$ in $C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)$. Now, we can see that $U_n-U$ satisfies

(11)$\left\{\begin{array}{cc}{\left(U_n-U\right)}_t=\Delta \left(U_n-U\right)+f_n,\hfill & x\in \Omega ,t\in \left(0,T_{max}\right),\hfill \\ {\displaystyle \frac{\partial (U_n-U)}{\partial n}}+g_1(u_n)u_n-g_1(u)u=0,\hfill \end{array}\right.$

where $f_n:=k(x,t){\int }_{\Omega }(f_1(u_n,v_n)-f_1(u,v))dx$. It is not difficult to verify that $f_n$ satisfies the assumptions of Theorem 2. We claim that $U_n\to U$ in $C^{2+\alpha }\left({\overline{\Omega }}_T\right)$ by using (8). Similarly, we can obtain that $V_n\to V$ in $C^{2+\alpha }\left({\overline{\Omega }}_T\right)$.

In order to apply the Leray–Schauder fixed-point theorem, we only need to prove that if T is sufficiently small, and $M\ge 2\left(1+d{(\Omega )}^{1-\alpha }\right)max\{{∥u_0∥}_{C^{2+\alpha }(\overline{\Omega })},{∥v_0∥}_{C^{2+\alpha }(\overline{\Omega })}\}$, then $A\left(B_T(M)\right)\subset B_T(M)$. A direct calculation shows that

$\begin{array}{c}\hfill \left|A_1(u(x,t),v(x,t))\right|\le \left|A_1(u(x,0),v(x,0))\right|+t{∥{\partial }_t{A}_1(u,v)∥}_{C^0(\overline{\Omega })}\le {∥u_0∥}_{C^0(\overline{\Omega })}+T{M}^{\prime },\end{array}$

$\begin{array}{cc}\hfill & \left|{\partial }_{x_i}{A}_1(u(x,t),v(x,t))\right|\hfill \\ \hfill \le & \left|{\partial }_{x_i}\left(A_1(u(x,0),v(x,0))\right)\right|+\mid {\partial }_{x_i}\left(A_1(u(x,t),v(x,t))-A_1(u(x,0),v(x,0))\mid \right.\hfill \\ \hfill \le & {∥{\partial }_{x_i}{u}_0∥}_{C^0(\overline{\Omega })}+t^{{\displaystyle \frac{1+\alpha }{2}}}{⟨A_1(u,v)⟩}_{2+\alpha ,{\overline{\Omega }}_T}\le {∥{\partial }_{x_i}{u}_0∥}_{C^0(\overline{\Omega })}+M^{\prime }{T}^{{\displaystyle \frac{1+\alpha }{2}}},\hfill \end{array}$

$\begin{array}{cc}\hfill & {\displaystyle \frac{\mid A_1(u(x,t),v(x,t))-A_1(u(x,s),v(x,s)\mid }{{|t-s|}^{\frac{1+\alpha }{2}}}}\le {∥{\partial }_t{A}_1(u,v)∥}_{C^0\left({\overline{\Omega }}_T\right)}{|t-s|}^{{\displaystyle \frac{1-\alpha }{2}}}\le M^{\prime }{T}^{{\displaystyle \frac{1-\alpha }{2}}},\hfill \end{array}$

and

$\begin{array}{cc}\hfill & \left|{\partial }_{x_i{x}_j}{A}_1(u(x,t),v(x,t))\right|\le \left|{\partial }_{x_i{x}_j}{u}_0(x)\right|+t^{{\displaystyle \frac{\alpha }{2}}}{\left[A_1(u,v)\right]}_{2+\alpha ,{\overline{\Omega }}_T}\le {∥{\partial }_{x_i{x}_j}{u}_0∥}_{C^0(\overline{\Omega })}+M^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}.\hfill \end{array}$

Combined with $\Omega$, which is a convex set, we can obtain

$\left|{\partial }_{x_i}{A}_1(u(x,s),v(x,s))-{\partial }_{x_i}{A}_1(u(y,s),v(y,s))\right|\le \parallel x-y\parallel \left(\sum _j{∥{\partial }_{x_i{x}_j}{u}_0∥}_{C^0(\overline{\Omega })}+M^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}\right).$

It follows that

$\begin{array}{cc}\hfill & \left|{\partial }_{x_i}{A}_1(u(x,t),v(x,t))-{\partial }_{x_i}{A}_1(u(y,s),v(y,s))\right|\hfill \\ \hfill \le & \left|{\partial }_{x_i}{A}_1(u(x,t),v(x,t))-{\partial }_{x_i}{A}_1(u(x,s),v(x,s))\right|\hfill \\ \hfill +& \left|{\partial }_{x_i}{A}_1(u(x,s),v(x,s))-{\partial }_{x_i}{A}_1(u(y,s),v(y,s))\right|\hfill \\ \hfill \le & M^{\prime }{|t-s|}^{{\displaystyle \frac{1+\alpha }{2}}}+\parallel x-y\parallel \left(\sum _j{∥{\partial }_{x_i{x}_j}{u}_0∥}_{C^0(\overline{\Omega })}+M^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}\right).\hfill \end{array}$

This implies

$\begin{array}{cc}\hfill {\displaystyle \frac{\left|{\partial }_{x_i}{A}_1(u(x,t),v(x,t))-{\partial }_{x_i}{A}_1(u(y,s),v(y,s))\right|}{{|x-y|}^{\alpha }+{|t-s|}^{\frac{\alpha }{2}}}}& \le d{(\Omega )}^{1-\alpha }\left(M^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}+\sum _j{∥{\partial }_{x_i{x}_j}{u}_0∥}_{C^0(\overline{\Omega })}\right)\hfill \\ \hfill & +M^{\prime }{T}^{{\displaystyle \frac{1}{2}}}.\hfill \end{array}$

Using these estimates, we conclude that

$\begin{array}{cc}\hfill \parallel A_1{(u,v)\parallel }_{C^{1+\alpha }\left({\overline{\Omega }}_T\right)}& \le \left(1+d{(\Omega )}^{1-\alpha }\right){\parallel u_0\parallel }_{C^0(\overline{\Omega })}+T{M}^{\prime }+n{T}^{{\displaystyle \frac{1+\alpha }{2}}}{M}^{\prime }\hfill \\ \hfill & +T^{{\displaystyle \frac{1-\alpha }{2}}}{M}^{\prime }+n{M}^{\prime }{T}^{{\displaystyle \frac{1}{2}}}+nd{(\Omega )}^{1-\alpha }{M}^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}.\hfill \end{array}$

Since the two components of the map $A(u,v):=(A_1(u,v),A_2(u,v))=(U,V)$ are symmetric, by repeating the above calculation, we can similarly obtain

$\begin{array}{cc}\hfill \parallel A_2{(u,v)\parallel }_{C^{1+\alpha }\left({\overline{\Omega }}_T\right)}& \le \left(1+d{(\Omega )}^{1-\alpha }\right){\parallel u_0\parallel }_{C^0(\overline{\Omega })}+T{M}^{\prime }+n{T}^{{\displaystyle \frac{1+\alpha }{2}}}{M}^{\prime }\hfill \\ \hfill & +T^{{\displaystyle \frac{1-\alpha }{2}}}{M}^{\prime }+n{M}^{\prime }{T}^{{\displaystyle \frac{1}{2}}}+nd{(\Omega )}^{1-\alpha }{M}^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}.\hfill \end{array}$

Then,

$\begin{array}{cc}\hfill {\parallel A(u,v)\parallel }_{C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)}& \le \left(1+d{(\Omega )}^{1-\alpha }\right)max\{\parallel u_0{\parallel }_{C^0(\overline{\Omega })},\parallel v_0{\parallel }_{C^0(\overline{\Omega })}\}\hfill \\ \hfill & +T{M}^{\prime }+n{T}^{{\displaystyle \frac{1+\alpha }{2}}}{M}^{\prime }+T^{{\displaystyle \frac{1-\alpha }{2}}}{M}^{\prime }+n{M}^{\prime }{T}^{{\displaystyle \frac{1}{2}}}+nd{(\Omega )}^{1-\alpha }{M}^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}.\hfill \end{array}$

Since $M^{\prime }$ is independent of T for all $T<1$, we can choose T that is sufficiently small such that

$T{M}^{\prime }+n{T}^{{\displaystyle \frac{1+\alpha }{2}}}{M}^{\prime }+T^{{\displaystyle \frac{1-\alpha }{2}}}{M}^{\prime }+n{M}^{\prime }{T}^{{\displaystyle \frac{1}{2}}}+nd{(\Omega )}^{1-\alpha }{M}^{\prime }{T}^{{\displaystyle \frac{\alpha }{2}}}\le {\displaystyle \frac{M}{2}}.$

This further implies that

${\parallel A(u,v)\parallel }_{C^{1+\alpha }\left({\overline{\Omega }}_T\right)\times C^{1+\alpha }\left({\overline{\Omega }}_T\right)}\le M\mathrm{for}\mathrm{all}u\in B_T(M).$

Thus, A has a fixed point in $B_T(M)$.

Now, if $(u,v)$ is a fixed point of $(A_1,A_2)$, then $(u,v)\in C^{2+\alpha }\left({\overline{\Omega }}_T\right)\times C^{2+\alpha }\left({\overline{\Omega }}_T\right)$ is a solution of (1).

The uniqueness of classical solutions will now be proved by a contradiction. Suppose that $\left(u_1,v_1\right)$ and $\left(u_2,v_2\right)$ are two classical solutions of (1). Let $E:=u_1-u_2$, $\widehat{E}:=v_1-v_2$. Then, $(E,\widehat{E})$ is a solution of the following system:

(12)$\left\{\begin{array}{cc}{E}_t=\Delta E+L_1\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ {\widehat{E}}_t=\Delta \widehat{E}+L_2\hfill & (x,t)\in \Omega \times (0,T],\hfill \end{array}\right.$

where

(13)$\begin{array}{cc}\hfill & L_1=k_1(x,t){\int }_{\Omega }{f}_1\left(u_1,v_1\right)-f_1\left(u_2,v_2\right)\mathrm{d}x,\hfill \\ \hfill & L_2=k_2(x,t){\int }_{\Omega }{f}_2\left(u_1,v_1\right)-f_2\left(u_2,v_2\right)\mathrm{d}x\hfill \end{array}$

and the boundary condition

(14)${\displaystyle \frac{\partial E}{\partial \mathit{n}}}=g_1(u_2)u_2-g_1(u_1)u_1,{\displaystyle \frac{\partial \widehat{E}}{\partial \mathit{n}}}=g_2(u_2)u_2-g_2(u_1)u_1,(x,t)\in \partial \Omega \times (0,T].$

Multiplying the first equation of (12) by E, and using Lemma 2 and the trace theorem, we obtain

(15)$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\parallel E\parallel }^2+{\parallel \nabla E\parallel }^2& \le C|E|\left|g_1\left(u_1\right)u_1-g_1\left(u_2\right)u_2\right|+C\parallel E\parallel ∥k_1{\int }_{\Omega }{f}_1\left(u_1,v_1\right)-f_1\left(u_2,v_2\right)\mathrm{d}x∥\hfill \\ \hfill & \le {C\epsilon \parallel \nabla E\parallel }^2+{(C+C(\epsilon ))\parallel E\parallel }^2+C{\parallel \widehat{E}\parallel }^2\hfill \end{array}$

${\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{∥E∥}^2+\left(1-C\epsilon \right){∥\nabla E∥}^2\le \left(C+C\left(\epsilon \right)\right){∥E∥}^2+C{\parallel \widehat{E}\parallel }^2.$

Take $\epsilon ={\epsilon }_0$ to be sufficiently small such that $1-C{\epsilon }_0>0$, and then integrate the time variable from 0 to t

(16)${∥E∥}^2\le {\int }_0^t({∥E∥}^2+\parallel \widehat{E}{\parallel }^2)\mathrm{d}s+{∥E\left(0\right)∥}^2+{\parallel \widehat{E}(0)\parallel }^2.$

Using a similar discussion, we can also obtain

(17)${∥\widehat{E}∥}^2\le {\int }_0^t({∥E∥}^2+\parallel \widehat{E}{\parallel }^2)\mathrm{d}s+{∥E\left(0\right)∥}^2+{\parallel \widehat{E}(0)\parallel }^2.$

Note that $\parallel E\left(0\right){\parallel }^2={\parallel \widehat{E}\left(0\right)\parallel }^2=0$. By adding (16) and (17), and applying Gronwall’s lemma, we obtain the uniqueness of the classical solution.

□

3. L² Estimate for Continuous-Time Galerkin Approximation

In this section, we propose the continuous-time Galerkin approximation scheme for Equation (1) and establish the uniqueness of the approximation solution. Subsequently, we obtain the optimal-order error estimation in $L^2$ norm by employing the nonlinear projection $\mathcal{U}(t)$ and $\mathcal{V}(t)$.

Let $\left\{S_h\right\}$ be a family of finite-dimensional subspaces of $H^1(\Omega )$ related to parameter h with the following properties:

For some $r\ge 1,\alpha =0,1$, and $\alpha \le \beta \le r+1$, there exists a constant $C>0$ such that

(18)$\underset{\chi \in S_h}{inf}{\parallel w-\chi \parallel }_{H^{\alpha }(\Omega )}\le C{\parallel w\parallel }_{H^{\beta }(\Omega )}{h}^{\beta -\alpha },w\in H^{\beta }(\Omega ).$

Moreover, we have

(19)${\parallel \xi \parallel }_{L^{\infty }(\Omega )}\le C{h}^{-d/2}{\parallel \xi \parallel }_{L^2(\Omega )}{,\parallel \xi \parallel }_{H^1(\Omega )}\le C{h}^{-1}{\parallel \xi \parallel }_{L^2(\Omega )},\mathrm{for}\xi \in S_h.$

Assumptions (18) and (19) are standard ones that are satisfied by the approximation spaces of piecewise polynomials defined over a regular family of triangulations.

To obtain the optimal error estimate for the Galerkin approximate solution of Equation (1), we further assume that $d\le 3$, $g_i\ge 0(i=1,2)$ and $g_i\in C^3(\Omega )$. In addition, both u and v satisfy the following regularity conditions:

(20)$\begin{array}{cc}\hfill & u,v\in L^{\infty }\left(0,T,H^{r+1}(\Omega )\right),\hfill \\ \hfill & u_t,v_t\in L^2\left(0,T,H^{r+1}(\Omega )\right)\cap L^{\infty }\left(0,T,L^{\infty }(\Omega )\right)\cap L^{\infty }(0,T,H^{{\displaystyle \frac{3}{2}}}(\partial \Omega )),\hfill \\ \hfill & \nabla u_t,\nabla v_t\in L^1\left(0,T,L^{\infty }(\Omega )\right),u_{tt},v_{tt}\in L^{\infty }\left(0,T,H^1(\Omega )\right),\mathrm{and}\hfill \\ \hfill & u_{ttt},v_{ttt}\in L^1\left(0,T,H^1(\Omega )\right)\cap L^2\left(0,T,L^2(\Omega )\right).\hfill \end{array}$

It can be seen that Condition (20) implies that $u_0,v_0$ have higher-order continuous derivatives.

Our analysis builds upon the nonlinear elliptic $H^1$ projection introduced in [12], which provides an effective framework for the boundary norm associated with nonlinear boundary conditions. Let $\mathcal{U}(t),\mathcal{V}(t):[0,T]\to S_h$ be the nonlinear $H^1$ projections of $u(t),v(t)$, respectively, defined by the following formula:

(21)$\begin{array}{cc}(\nabla (u(t)-\mathcal{U}(t)),\nabla \varphi )+\lambda (u(t)-\mathcal{U}(t),\varphi )+⟨g_1(u(t))u(t)-g_1(\mathcal{U}(t))\mathcal{U}(t),\varphi ⟩=0,\hfill & \varphi \in S_h,\hfill \\ (\nabla (v(t)-\mathcal{V}(t)),\nabla \psi )+\lambda (v(t)-\mathcal{V}(t),\psi )+⟨g_2(v(t))v(t)-g_2(\mathcal{V}(t))\mathcal{V}(t),\psi ⟩=0,\hfill & \psi \in S_h\hfill \end{array}$

where $\lambda$ is a sufficiently large positive constant that guarantees the existence and uniqueness of $\mathcal{U}(t)$ and $\mathcal{V}(t)$. The approximation error between the exact solution of the equation and the nonlinear projection is described by the following result:

Lemma 3 

(see [12]). Let $\mathcal{U}(t)$ and $\mathcal{V}(t)$ be the projections defined by the above formula. We set $\rho =u-\mathcal{U},\widehat{\rho }=v-\mathcal{V}$. Then, there exists a positive constant C such that the following inequality holds:

(22)$\begin{array}{c}{\parallel \rho \parallel }_{L^{\infty }\left(L^2\right)}+\parallel \widehat{\rho }{\parallel }_{L^{\infty }\left(L^2\right)}+{h\parallel \rho \parallel }_{L^{\infty }\left(H^1\right)}+h{\parallel \widehat{\rho }\parallel }_{L^{\infty }\left(H^1\right)}\le C_1{h}^{r+1},\\ {∥{\rho }_t∥}_{L^2\left(L^2\right)}+{∥{\widehat{\rho }}_t∥}_{L^2\left(L^2\right)}\le C_1{h}^{r+1},\\ {∥{\rho }_t∥}_{L^{\infty }\left(H^1\right)}+{∥{\rho }_{tt}∥}_{L^{\infty }\left(H^1\right)}+{∥{\rho }_{ttt}∥}_{L^1\left(H^1\right)}\le C_1,\\ {∥{\widehat{\rho }}_t∥}_{L^{\infty }\left(H^1\right)}+{∥{\widehat{\rho }}_{tt}∥}_{L^{\infty }\left(H^1\right)}+{∥{\widehat{\rho }}_{ttt}∥}_{L^1\left(H^1\right)}\le C_1.\end{array}$

In fact, $\rho$ and $\widehat{\rho }$ have exactly the same properties, because the projections $\mathcal{U}(t)$ and $\mathcal{V}(t)$ are defined in the same way. We also assume that the constant $C_1$ can be selected to satisfy the following requirements:

(23)$\begin{array}{cc}\hfill & {\parallel \mathcal{U}\parallel }_{L^{\infty }\left(L^{\infty }\right)}+{\parallel \mathcal{V}\parallel }_{L^{\infty }\left(L^{\infty }\right)}+{\parallel \nabla \mathcal{U}\parallel }_{L^{\infty }\left(L^{\infty }\right)}+{\parallel \nabla \mathcal{V}\parallel }_{L^{\infty }\left(L^{\infty }\right)}\le C_1,\hfill \\ \hfill & {∥{\mathcal{U}}_t∥}_{L^1\left(L^{\infty }\right)}+{∥{\mathcal{V}}_t∥}_{L^1\left(L^{\infty }\right)}+{∥\nabla {\mathcal{U}}_t∥}_{L^1\left(L^{\infty }\right)}+{∥\nabla {\mathcal{V}}_t∥}_{L^1\left(L^{\infty }\right)}\le C_1.\hfill \end{array}$

In fact, this assumption is valid only in the case of $d=3,r=1$. For $d=3,r\ge 2$ or $d\le 2$, it can be derived from Equations (18), (19) and (22).

The continuous-time Galerkin approximation of System (1) is to find $\left\{U\left(t\right),V\left(t\right)\right\}:\left[0,T\right]\to S_h\times S_h$, such that

(24)$\left\{\begin{array}{cc}\left(U_t,\varphi \right)+(\nabla U,\nabla \varphi )+⟨g_1(U)U,\varphi ⟩=\left(k_1{\int }_{\Omega }{f}_1\left(U,V\right)\mathrm{d}x,\varphi \right),\hfill & \forall \varphi \in S_h,\hfill \\ \left(V_t,\psi \right)+(\nabla V,\nabla \psi )+⟨g_2(V)V,\psi ⟩=\left(k_2{\int }_{\Omega }{f}_2\left(U,V\right)\mathrm{d}x,\psi \right),\hfill & \forall \psi \in S_h.\hfill \end{array}\right.$

First, we consider the uniqueness of approximation $U(t),V(t)$.

Theorem 3. 

The approximations $U(t),V(t)$ defined by Equation (24) are unique.

Proof. 

The theorem can be derived by using the proof method of the uniqueness of the classical solution in Section 2. □

For the error between the true solution and the approximate solution, we have

Theorem 4. 

Let $U\left(t\right)$ and $V\left(t\right)$ be continuous Galerkin approximation solutions satisfying $∥U\left(·,0\right)-u_0∥+∥V\left(·,0\right)-v_0∥\le C{h}^{r+1}$. Thus, we have

$\begin{array}{c}\hfill {∥u-U∥}_{L^{\infty }\left(L^2\right)}+{∥v-V∥}_{L^{\infty }\left(L^2\right)}\le C{h}^{r+1}.\end{array}$

Proof. 

First, we decompose the error as follows:

$u-U=(u-\mathcal{U})+(\mathcal{U}-U)=\rho +\sigma ,v-V=(v-\mathcal{V})+(\mathcal{V}-V)=\widehat{\rho }+\widehat{\sigma }.$

Combining Equations (1) and (24), we can obtain

(25)$\begin{array}{cc}\hfill & \left({\mathcal{U}}_t,\varphi \right)+(\nabla \mathcal{U},\nabla \varphi )+⟨g_1(\mathcal{U})\mathcal{U},\varphi ⟩=\left(k_1{\int }_{\Omega }{f}_1(u,v)\mathrm{d}x-{\rho }_t,\varphi \right)+\lambda (\rho ,\varphi ),\hfill \\ \hfill & \left({\mathcal{V}}_t,\psi \right)+(\nabla \mathcal{V},\nabla \psi )+⟨g_2(\mathcal{V})\mathcal{V},\psi ⟩=\left(k_2{\int }_{\Omega }{f}_2(u,v)\mathrm{d}x-{\widehat{\rho }}_t,\psi \right)+\lambda (\widehat{\rho },\psi ).\hfill \end{array}$

Subtracting Equation (25) with (24), we obtain

(26)$\begin{array}{cc}\hfill & \left({\sigma }_t,\varphi \right)+(\nabla \sigma ,\nabla \varphi )+⟨g_1(\mathcal{U})\mathcal{U}-g_1(U)U,\varphi ⟩\hfill \\ \hfill =& \left(k_1{\int }_{\Omega }\left(f_1(u,v)-f_1(U,V)\right)\mathrm{d}x-{\rho }_t,\varphi \right)-\lambda (\rho ,\varphi ),\hfill \end{array}$

(27)$\begin{array}{cc}\hfill & \left({\widehat{\sigma }}_t,\psi \right)+(\nabla \widehat{\sigma },\nabla \psi )+⟨g_2(\mathcal{V})\mathcal{V}-g_2(V)V,\psi ⟩\hfill \\ \hfill =& \left(k_2{\int }_{\Omega }\left(f_2(u,v)-f_2(U,V)\right)\mathrm{d}x-{\widehat{\rho }}_t,\psi \right)-\lambda (\widehat{\rho },\psi ).\hfill \end{array}$

Taking $\varphi =\sigma$ in (26), we obtain

$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\parallel \sigma \parallel }^2+{\parallel \nabla \sigma \parallel }^2-C{|\sigma |}^2\hfill \\ \hfill & \le {C(\parallel \rho \parallel }^2+{\parallel \sigma \parallel }^2+{∥{\rho }_t∥}^2)+C{∥k_1{\int }_{\Omega }(f_1(u,v)-f_1(U,v)+f_1(U,v)-f_1(U,V))\mathrm{d}x∥}^2.\hfill \end{array}$

By using Lemma 2 and the trace theorem, we have

$\begin{array}{c}\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\parallel \sigma \parallel }^2+{\parallel \nabla \sigma \parallel }^2\le {C(\parallel \rho \parallel }^2+{\parallel \sigma \parallel }^2+{∥{\rho }_t∥}^2+\parallel \widehat{\rho }{\parallel }^2+\parallel \widehat{\sigma }{\parallel }^2)+{\displaystyle \frac{C_1}{8}}{\parallel \nabla \sigma \parallel }^2.\hfill \end{array}\hfill \end{array}$

According to Hypothesis (22), there is

$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\parallel \sigma \parallel }^2+{\displaystyle \frac{C_1}{2}}{\parallel \nabla \sigma \parallel }^2\le C\left({\parallel \sigma \parallel }^2+{\parallel \widehat{\sigma }\parallel }^2\right)+C{h}^{2r+2}.\hfill \end{array}$

The following result can be obtained by integrating both sides of the above inequality from 0 to t:

${\parallel \sigma \parallel }^2-{\parallel \sigma (0)\parallel }^2+C_1{\int }_0^t{\parallel \nabla \sigma \parallel }^2\mathrm{d}s\le C{\int }_0^t\left({\parallel \sigma \parallel }^2+{\parallel \widehat{\sigma }\parallel }^2\right)\mathrm{d}s+C{h}^{2r+2}.$

Taking $\psi =\widehat{\sigma }$ in (27), and similar to the above discussion, we can also obtain

$\parallel \widehat{\sigma }{\parallel }^2-\parallel \widehat{\sigma }{(0)\parallel }^2+C_1{\int }_0^t{\parallel \nabla \widehat{\sigma }\parallel }^2\mathrm{d}s\le C{\int }_0^t\left({\parallel \sigma \parallel }^2+{\parallel \widehat{\sigma }\parallel }^2\right)\mathrm{d}s+C{h}^{2r+2}.$

Thus,

(28)$\begin{array}{c}\begin{array}{cc}\hfill & {\parallel \sigma \parallel }^2+\parallel \widehat{\sigma }{\parallel }^2\le C{\int }_0^t\left({\parallel \sigma \parallel }^2+{\parallel \widehat{\sigma }\parallel }^2\right)\mathrm{d}s+C{h}^{2r+2}+{\parallel \sigma (0)\parallel }^2+{\parallel \widehat{\sigma }(0)\parallel }^2.\hfill \end{array}\hfill \end{array}$

Notice that

$\begin{array}{cc}\hfill & {\parallel \sigma (0)\parallel }^2\le C({∥u_0-U(·,0)∥}^2+{\parallel \rho (0)\parallel }^2)\le C{h}^{2r+2},\hfill \\ \hfill & {\parallel \sigma (0)\parallel }^2\le C({∥u_0-U(·,0)∥}^2+{\parallel \rho (0)\parallel }^2)\le C{h}^{2r+2}.\hfill \end{array}$

Then, by (28), Gronwall’s lemma, and the triangle inequality, the conclusion holds. □

4. L² Estimate for Crank–Nicolson Galerkin Approximation

In this section, we discretize the time interval and provide the Crank–Nicolson Galerkin approximation of System (1). The uniqueness of the approximation solution and the optimal-order error estimate in $L^2$ norm are obtained.

Before presenting the main results of this section, we provide some necessary explanations for some notations. Let N be a positive integer, and let $0=t_0\le t_1\le t_2\le \cdots \le t_N=T$ be a partition of the interval $[0,T]$, where the time step $▵t={\displaystyle \frac{T}{N}}$. Assume that $\varphi (x,t)$ is a function defined on $\Omega \times {\bigcup }_{n=1}^N\left\{t_n\right\}$. Thus, we define

$\begin{array}{cc}\hfill & t_n=n\Delta t,t_{n+{\displaystyle \frac{1}{2}}}=\left(n+{\displaystyle \frac{1}{2}}\right)\Delta t,{\varphi }^n=\varphi \left(x,t_n\right),\hfill \\ \hfill & \partial {\varphi }^n={\varphi }^{n+1}-{\varphi }^n,{\varphi }^{n+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{{\varphi }^{n+1}+{\varphi }^n}{2}},{\partial }_t{\varphi }^n={\displaystyle \frac{{\varphi }^{n+1}-{\varphi }^n}{\Delta t}},\hfill \end{array}$

and we specifically note that, in general, ${\varphi }^{n+{\displaystyle \frac{1}{2}}}\ne \varphi (x,t_{n+{\displaystyle \frac{1}{2}}})$.

The Crank–Nicolson Galerkin approximation of System (1) is to find sequences ${\left\{U_n,V_n\right\}}_{n=0}^N\in S_h\times S_h$, such that

(29)$\left\{\begin{array}{cc}\left({\partial }_t{U}_n,\varphi \right)\hfill & +(\nabla U_{n+{\displaystyle \frac{1}{2}}},\nabla \varphi )+⟨g_1(U_{n+{\displaystyle \frac{1}{2}}})U_{n+{\displaystyle \frac{1}{2}}},\varphi ⟩\hfill \\ & =(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x,\varphi )\forall \varphi \in S_h,\hfill \\ \left({\partial }_t{V}_n,\psi \right)\hfill & +(\nabla V_{n+{\displaystyle \frac{1}{2}}},\nabla \psi )+⟨g_2(V_{n+{\displaystyle \frac{1}{2}}})V_{n+{\displaystyle \frac{1}{2}}},\psi ⟩\hfill \\ & =(k_2(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_2(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x,\psi )\forall \psi \in S_h.\hfill \end{array}\right.$

Theorem 5. 

The sequences ${\left\{U_n,V_n\right\}}_{n=0}^N$ defined by Equation (29) exist and are unique.

Proof. 

First, we prove the existence of fully discrete solutions.

We proceed by induction, assuming that the sequences ${\left\{U_k,V_k\right\}}_{k=0}^n$ are given. For $\mathit{z}=\left(z_1,z_2\right)\in S_h\times S_h$, we define a mapping $F:\left(S_h\times S_h\right)\to \left(S_h\times S_h\right)$ such that

(30)$\begin{array}{cc}\hfill \left(\mathit{F}\left(\mathit{z}\right),\Psi \right)& =\left(\mathit{z}-{\mathit{x}}_{\mathit{n}},\Psi \right)+{\displaystyle \frac{1}{2}}\Delta t\left(\nabla \mathit{z},\nabla \Psi \right)+{\displaystyle \frac{1}{2}}\Delta t⟨\mathit{g}(\mathit{z})\mathit{z},\Psi ⟩\hfill \\ & -{\displaystyle \frac{1}{2}}\Delta t\left(\mathit{k}{\int }_{\Omega }\mathit{f}\left(z_1,z_2\right)\mathrm{d}x,\Psi \right),\forall \Psi \in S_h\times S_h\hfill \end{array}$

where $\Psi =\left(\varphi ,\psi \right),{\mathit{x}}_{\mathit{n}}=\left(U_n,V_n\right),\mathit{g}\left(\mathit{z}\right)\mathit{z}=\left(g_1\left(z_1\right)z_1,g_2\left(z_2\right)z_2\right),\mathit{f}=\left(f_1\left(z_1,z_2\right),f_2\left(z_1,z_2\right)\right),\mathit{k}=\left(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}),k_2(x,t_{n+{\displaystyle \frac{1}{2}}})\right)$.

Taking $\Psi =\mathit{z}$ in the above equation and using the non-negativity of $g_i(i=1,2)$, we obtain

(31)$\begin{array}{cc}\hfill (F(\mathit{z}),\mathit{z})& =\left(\mathit{z}-{\mathit{x}}_{\mathit{n}},\mathit{z}\right)+{\displaystyle \frac{1}{2}}\Delta t{\parallel \nabla \mathit{z}\parallel }^2+{\displaystyle \frac{1}{2}}\Delta t⟨\mathit{g}(\mathit{z})\mathit{z},\mathit{z}⟩-{\displaystyle \frac{1}{2}}\Delta t\left(\mathit{k}{\int }_{\Omega }\mathit{f}\left(z_1,z_2\right)\mathrm{d}x,\mathit{z}\right)\hfill \\ \hfill & \ge {\parallel \mathit{z}\parallel }^2-∥{\mathit{x}}_{\mathit{n}}∥\parallel \mathit{z}\parallel -{\displaystyle \frac{1}{2}}\Delta t\parallel \mathit{k}(x,t_{n+{\displaystyle \frac{1}{2}}})\parallel {\int }_{\Omega }\mathit{f}\left(z_1,z_2\right)\mathrm{d}x\parallel \parallel \mathit{z}\parallel \hfill \\ \hfill & \ge {\parallel \mathit{z}\parallel }^2-∥{\mathit{x}}_{\mathit{n}}∥\parallel \mathit{z}\parallel -{\displaystyle \frac{1}{2}}C\Delta t∥\mathit{f}\left(z_1,z_2\right)∥\parallel \mathit{z}\parallel \hfill \\ \hfill & \ge {\displaystyle \frac{1}{2}}\left({\parallel \mathit{z}\parallel }^2-{∥{\mathit{x}}_{\mathit{n}}∥}^2\right)-{\displaystyle \frac{1}{2}}C\Delta t{\parallel \mathit{z}\parallel }^2.\hfill \end{array}$

Taking $\Delta t$ to be sufficiently small such that ${\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{2}}C\Delta t\ge {\displaystyle \frac{1}{4}}$, and choosing ${∥\mathit{z}∥}^2>2{∥{\mathit{x}}_{\mathit{n}}∥}^2$, we derive

$\left(F\left(\mathit{z}\right),\mathit{z}\right)\ge {\displaystyle \frac{1}{4}}\left({∥\mathit{z}∥}^2-2{∥{\mathit{x}}_{\mathit{n}}∥}^2\right)>0.$

By Brouwer’s fixed-point theorem, there exists ${\mathit{z}}^{*}=\left(z_1^{*},z_2^{*}\right)\in S_h\times S_h$ such that $F\left({\mathit{z}}^{*}\right)=0$. Let $U_{n+1}=2z_1^{*}-U_n$ and $V_{n+1}=2z_2^{*}-V_n$. It is easy to verify that $U_{n+1},V_{n+1}$ is a solution satisfying (29). Therefore, we conclude the existence of the sequence of solutions.

Next, we consider the uniqueness of the solution sequence. First, choose $U_0={\widehat{U}}_0,V_0={\widehat{V}}_0$ such that $E_0={\widehat{E}}_0=0$. Moreover, we only need to prove that $E_{n+1}={\widehat{E}}_{n+1}=0$ under the condition $E_n={\widehat{E}}_n=0$.

Using the same notation as in Theorem 3, and following a similar method, we have

(32)$\begin{array}{cc}\hfill ({\partial }_t{E}_n,\varphi )& +(\nabla E_{n+{\displaystyle \frac{1}{2}}},\nabla \varphi )+⟨g_1(U_{n+{\displaystyle \frac{1}{2}}})U_{n+{\displaystyle \frac{1}{2}}}-g_1({\widehat{U}}_{n+{\displaystyle \frac{1}{2}}}){\widehat{U}}_{n+{\displaystyle \frac{1}{2}}},\varphi ⟩\hfill \\ \hfill & =\left(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})-f_1({\widehat{U}}_{n+{\displaystyle \frac{1}{2}}},{\widehat{V}}_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x,\varphi \right),\hfill \end{array}$

(33)$\begin{array}{cc}\hfill ({\partial }_t{\widehat{E}}_n,\psi )& +(\nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}},\nabla \psi )+⟨g_2(V_{n+{\displaystyle \frac{1}{2}}})V_{n+{\displaystyle \frac{1}{2}}}-g_2({\widehat{V}}_{n+{\displaystyle \frac{1}{2}}}){\widehat{V}}_{n+{\displaystyle \frac{1}{2}}},\psi ⟩\hfill \\ \hfill & =\left(k_2(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_2(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})-f_2({\widehat{U}}_{n+{\displaystyle \frac{1}{2}}},{\widehat{V}}_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x,\psi \right).\hfill \end{array}$

Taking $\varphi =E_{n+{\displaystyle \frac{1}{2}}}$ in (32), by Lemma 2 and the trace theorem, we obtain

$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2\Delta t}}(\parallel E_{n+1}{\parallel }^2-\parallel E_n{\parallel }^2)+\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2\hfill \\ \hfill \le & C|g_1(U_{n+{\displaystyle \frac{1}{2}}})U_{n+{\displaystyle \frac{1}{2}}}-g_1({\widehat{U}}_{n+{\displaystyle \frac{1}{2}}}){\widehat{U}}_{n+{\displaystyle \frac{1}{2}}}{|}^2+C{|E_{n+{\displaystyle \frac{1}{2}}}|}^2\hfill \\ \hfill +& C{∥{\int }_{\Omega }{f}_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})-f_1({\widehat{U}}_{n+{\displaystyle \frac{1}{2}}},{\widehat{V}}_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x∥}^2+C{\parallel E_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill \le & C|E_{n+{\displaystyle \frac{1}{2}}}{|}^2+C\parallel f_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})-f_1({\widehat{U}}_{n+{\displaystyle \frac{1}{2}}},{\widehat{V}}_{n+{\displaystyle \frac{1}{2}}}){\parallel }^2+C{\parallel E_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill \le & C\epsilon \parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+(C(\epsilon )+C)\parallel E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+C{\parallel {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}\parallel }^2.\hfill \end{array}$

Then, multiplying by $2\Delta t$, we obtain

$\begin{array}{cc}\hfill & \parallel E_{n+1}{\parallel }^2-\parallel E_n{\parallel }^2+2\Delta t{\parallel E_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill \le & 2C\Delta t\epsilon \parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+2\Delta t(C(\epsilon )+C)\parallel E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+2C\Delta t{\parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}\parallel }^2.\hfill \end{array}$

Taking $\psi ={\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}$ in (33) and using a similar discussion, we obtain

$\begin{array}{cc}\hfill & \parallel {\widehat{E}}_{n+1}{\parallel }^2-\parallel {\widehat{E}}_n{\parallel }^2+2\Delta t{\parallel {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill \le & 2C\Delta t\epsilon \parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+2\Delta t(C(\epsilon )+C)\parallel {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+2C\Delta t{\parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}\parallel }^2.\hfill \end{array}$

Adding the two equations above and using the inequality $\parallel E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2\le {\displaystyle \frac{1}{2}}(\parallel E_n{\parallel }^2+\parallel E_{n+1}{\parallel }^2)$, we have

$\begin{array}{c}\begin{array}{cc}\hfill & (\parallel E_{n+1}{\parallel }^2+\parallel {\widehat{E}}_{n+1}{\parallel }^2)-(\parallel E_n\hfill \end{array}{\parallel }^2+\parallel {\widehat{E}}_n{\parallel }^2)+2\Delta t(\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \\ \hspace{1em}\hspace{1em}\le 2C\Delta t\epsilon (\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)+(C(\epsilon )+C)\Delta t(\parallel E_n{\parallel }^2+\parallel E_{n+1}{\parallel }^2+\parallel {\widehat{E}}_n{\parallel }^2+\parallel {\widehat{E}}_{n+1}{\parallel }^2),\hfill \end{array}$

thus,

$\begin{array}{cc}\hfill (1-\Delta t(C(\epsilon )+C))(\parallel E_{n+1}{\parallel }^2+\parallel {\widehat{E}}_{n+1}{\parallel }^2)+& 2\Delta t(1-C\epsilon )(\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \\ \hfill \le & (1+(C(\epsilon )+C)\Delta t)(\parallel E_n{\parallel }^2+\parallel {\widehat{E}}_n{\parallel }^2).\hfill \end{array}$

By using the arbitrariness of $\epsilon$, first take ${\epsilon }_0$ to be sufficiently small such that $2\Delta t(1-C{\epsilon }_0)>0$, and then fix ${\epsilon }_0$ and take $\Delta t$ to be sufficiently small such that $(1-(C({\epsilon }_0)+C)\Delta t)>0$. Then, there is

$\begin{array}{c}\begin{array}{cc}\hfill & \parallel E_{n+1}{\parallel }^2+\parallel {\widehat{E}}_{n+1}{\parallel }^2\le M(\parallel E_n{\parallel }^2+\parallel {\widehat{E}}_n{\parallel }^2)\hfill \end{array}\hfill \end{array}$

where $M={\displaystyle \frac{1+\Delta t(C+C({\epsilon }_0))}{1-\Delta t(C+C({\epsilon }_0))}}>0$. According to the inductive hypothesis, the conclusion is established. □

Theorem 6. 

Let ${\{U_n,V_n\}}_{n=0}^N$ be the approximate solution defined by (29). Then, for a sufficiently small $\Delta t$, we have

$\underset{0\le n\le N}{max}\parallel u_n-U_n\parallel +\parallel v_n-V_n\parallel \le C({(\Delta t)}^2+h^{r+1}),$

where C is a positive constant independent of h and $\Delta t$.

Proof. 

The notations used in the proof of Theorem 4 are still valid. By substituting $t=t_{n+{\displaystyle \frac{1}{2}}}$ in (21), we obtain the following for any $\varphi ,\psi \in S_h$:

(34)$\begin{array}{cc}\hfill & (\nabla u(t_{n+{\displaystyle \frac{1}{2}}}),\nabla \varphi )+⟨g_1(u(t_{n+{\displaystyle \frac{1}{2}}}))u(t_{n+{\displaystyle \frac{1}{2}}}),\varphi ⟩\hfill \\ \hfill =& (\nabla \mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}),\nabla \varphi )+⟨g_1(\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}))\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}),\varphi ⟩-\lambda (\rho (t_{n+{\displaystyle \frac{1}{2}}}),\varphi )\hfill \\ \hfill =& (\nabla {\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}},\nabla \varphi )+⟨g_1(\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}})\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}),\varphi ⟩-\lambda (\rho (t_{n+{\displaystyle \frac{1}{2}}}),\varphi )\hfill \\ \hfill +& (\nabla (\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}})-{\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}),\nabla \varphi ),\hfill \end{array}$

(35)$\begin{array}{cc}\hfill & (\nabla v(t_{n+{\displaystyle \frac{1}{2}}}),\nabla \psi )+⟨g_2(v(t_{n+{\displaystyle \frac{1}{2}}}))v(t_{n+{\displaystyle \frac{1}{2}}}),\psi ⟩\hfill \\ \hfill =& (\nabla \mathcal{V}(t_{n+{\displaystyle \frac{1}{2}}}),\nabla \psi )+⟨g_2(\mathcal{V}(t_{n+{\displaystyle \frac{1}{2}}}))\mathcal{V}(t_{n+{\displaystyle \frac{1}{2}}}),\psi ⟩-\lambda (\widehat{\rho }(t_{n+{\displaystyle \frac{1}{2}}}),\psi )\hfill \\ \hfill =& (\nabla {\mathcal{V}}_{n+{\displaystyle \frac{1}{2}}},\nabla \psi )+⟨g_2(\mathcal{V}(t_{n+{\displaystyle \frac{1}{2}}}))\mathcal{V}(t_{n+{\displaystyle \frac{1}{2}}}),\psi ⟩-\lambda (\widehat{\rho }(t_{n+{\displaystyle \frac{1}{2}}}),\psi )\hfill \\ \hfill +& (\nabla (\mathcal{V}(t_{n+{\displaystyle \frac{1}{2}}})-\nabla {\mathcal{V}}_{n+{\displaystyle \frac{1}{2}}}),\nabla \psi ).\hfill \end{array}$

We write the weak form of Equation (1) when $t=t_{n+{\displaystyle \frac{1}{2}}}$:

(36)$\begin{array}{cc}\hfill & (u_t(t_{n+{\displaystyle \frac{1}{2}}}),\varphi )+(\nabla u(t_{n+{\displaystyle \frac{1}{2}}}),\nabla \varphi )+⟨g_1(u(t_{n+{\displaystyle \frac{1}{2}}}))u(t_{n+{\displaystyle \frac{1}{2}}}),\varphi ⟩\hfill \\ \hfill =& \left(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_1(u(t_{n+{\displaystyle \frac{1}{2}}}),v(t_{n+{\displaystyle \frac{1}{2}}}))\mathrm{d}x,\varphi \right),\hfill \end{array}$

(37)$\begin{array}{cc}\hfill & (v_t(t_{n+{\displaystyle \frac{1}{2}}}),\psi )+(\nabla v(t_{n+{\displaystyle \frac{1}{2}}}),\nabla \psi )+⟨g_2(v(t_{n+{\displaystyle \frac{1}{2}}}))v(t_{n+{\displaystyle \frac{1}{2}}}),\psi ⟩\hfill \\ \hfill =& \left(k_2(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_2(u(t_{n+{\displaystyle \frac{1}{2}}}),v(t_{n+{\displaystyle \frac{1}{2}}}))\mathrm{d}x,\psi \right).\hfill \end{array}$

Combining (29), (34) and (36), we have

$\begin{array}{cc}\hfill & ({\partial }_t{\sigma }_n,\varphi )+(\nabla {\sigma }_{n+{\displaystyle \frac{1}{2}}},\nabla \varphi )\hfill \\ \hfill =& \left(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_1(u(t_{n+{\displaystyle \frac{1}{2}}}),v(t_{n+{\displaystyle \frac{1}{2}}}))-f_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x+{\eta }_n-{\partial }_t{\rho }_n,\varphi \right)\hfill \\ \hfill +& \lambda (\rho (t_{n+{\displaystyle \frac{1}{2}}}),\varphi )+(\nabla ({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}-\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}})),\nabla \varphi )+⟨g_1(U_{n+{\displaystyle \frac{1}{2}}})U_{n+{\displaystyle \frac{1}{2}}}-g_1({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}){\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}},\varphi ⟩\hfill \\ \hfill +& ⟨g_1({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}){\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}-g_1(\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}))\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}),\varphi ⟩\hfill \end{array}$

where ${\eta }_n={\partial }_t{u}_n-u_t(t_{n+{\displaystyle \frac{1}{2}}})$.

Setting $\varphi ={\sigma }_{k+{\displaystyle \frac{1}{2}}}\Delta t$, and using the Cauchy–Schwartz inequality, $\epsilon$—inequality and the trace theorem, we derive

(38)$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}& (\Vert {\sigma }_{n+1}{\Vert }^2-\Vert {\sigma }_n{\Vert }^2)+\Delta t\Vert \nabla {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\Vert }^2\hfill \\ \hfill \le & C\Delta t\Vert {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\Vert }^2+{C\Delta t\Vert \rho \Vert }_{L^{\infty }(L^2)}^2+C\epsilon \Delta t\Vert \nabla {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\Vert }^2+C(\epsilon )\Delta t{\Vert \nabla ({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}-\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}))\Vert }^2\hfill \\ \hfill +& C\Delta t\left({∥{\int }_{\Omega }{f}_1(u(t_{n+{\displaystyle \frac{1}{2}}}),v(t_{n+{\displaystyle \frac{1}{2}}})-f_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x∥}^2+\Vert {\eta }_n{\Vert }^2+{\Vert {\partial }_t{\rho }_n\Vert }^2\right)\hfill \\ \hfill +& C\Delta t|{\sigma }_{n+{\displaystyle \frac{1}{2}}}{|}^2+C\Delta t{|\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}})-{\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}|}^2\hfill \\ \hfill \le & (C+C(\epsilon ))\Delta t\Vert {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\Vert }^2+{C\Delta t\Vert \rho \Vert }_{L^{\infty }(L^2)}^2+C\epsilon \Delta t{\Vert \nabla {\sigma }_{n+{\displaystyle \frac{1}{2}}}\Vert }^2\hfill \\ \hfill +& C\Delta t\left({∥{\int }_{\Omega }{f}_1(u(t_{n+{\displaystyle \frac{1}{2}}}),v(t_{n+{\displaystyle \frac{1}{2}}}))-f_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x∥}^2+\Vert {\eta }_n{\Vert }^2+{\Vert {\partial }_t{\rho }_n\Vert }^2\right)\hfill \\ \hfill +& C(\epsilon )\Delta t(\Vert \nabla ({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}-\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}})){\Vert }^2+\Vert ({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}-\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}})){\Vert }^2).\hfill \end{array}$

Since

(39)$\begin{array}{cc}\hfill & {∥{\int }_{\Omega }{f}_1(u(t_{n+{\displaystyle \frac{1}{2}}}),v(t_{n+{\displaystyle \frac{1}{2}}}))-f_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x∥}^2\hfill \\ \hfill =& {∥{\int }_{\Omega }{f}_1(u(t_{n+{\displaystyle \frac{1}{2}}}),v(t_{n+{\displaystyle \frac{1}{2}}}))-f_1(U_{n+{\displaystyle \frac{1}{2}}},v(t_{n+{\displaystyle \frac{1}{2}}}))+f_1(U_{n+{\displaystyle \frac{1}{2}}},v(t_{n+{\displaystyle \frac{1}{2}}}))-f_1(U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}})\mathrm{d}x∥}^2\hfill \\ \hfill \le & C\parallel u(t_{n+{\displaystyle \frac{1}{2}}})-U_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+C{\parallel v(t_{n+{\displaystyle \frac{1}{2}}})-V_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill \le & C(\parallel u(t_{n+{\displaystyle \frac{1}{2}}})-u_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel u_{n+{\displaystyle \frac{1}{2}}}-U_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel v(t_{n+{\displaystyle \frac{1}{2}}})-v_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel v_{n+{\displaystyle \frac{1}{2}}}-V_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \\ \hfill \le & C(\parallel u(t_{n+{\displaystyle \frac{1}{2}}})-u_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel v(t_{n+{\displaystyle \frac{1}{2}}})-v_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \\ \hfill +& C(\parallel {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel {\rho }_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel {\widehat{\sigma }}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel {\widehat{\rho }}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \end{array}$

and

(40)$\Delta t\parallel {\partial }_t{\rho }_n{\parallel }^2\le {\int }_{t_n}^{t_{n+1}}{\parallel {\rho }_t\parallel }^2\mathrm{d}s.$

Using the Taylor formula with an integral remainder, the following inequality holds:

(41)$\begin{array}{cc}\hfill \parallel u(t_{n+{\displaystyle \frac{1}{2}}})-u_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+{\parallel v(t_{n+{\displaystyle \frac{1}{2}}})-v_{n+{\displaystyle \frac{1}{2}}}\parallel }^2& \le C{(\Delta t)}^3{\int }_{t_n}^{t_{n+1}}(\parallel u_{tt}{\parallel }^2+\parallel v_{tt}{\parallel }^2)\mathrm{d}s,\hfill \\ \hfill \parallel \nabla ({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}-\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}})){\parallel }^2+{\parallel ({\mathcal{U}}_{n+{\displaystyle \frac{1}{2}}}-\mathcal{U}(t_{n+{\displaystyle \frac{1}{2}}}))\parallel }^2& \le C{(\Delta t)}^3{\int }_{t_n}^{t_{n+1}}(\parallel \nabla {\mathcal{U}}_{tt}{\parallel }^2+\parallel {\mathcal{U}}_{tt}{\parallel }^2)\mathrm{d}s,\hfill \\ \hfill \parallel {\eta }_n{\parallel }^2& \le C{(\Delta t)}^3{\int }_{t_n}^{t_{n+1}}{\parallel u_{ttt}\parallel }^2\mathrm{d}s.\hfill \end{array}$

The following conclusion can be inferred from (38)–(41):

$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}(\parallel {\sigma }_{n+1}{\parallel }^2-\parallel {\sigma }_n{\parallel }^2)+(1-C\epsilon )\Delta t\parallel \nabla {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2\hfill \\ \hfill & \le C({(\Delta t)}^4+h^{2r+2})+(C+C(\epsilon ))\Delta t\parallel {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+C\Delta t{\parallel {\widehat{\sigma }}_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill & \le C({(\Delta t)}^4+h^{2r+2})+C(\epsilon )\Delta t\parallel {\sigma }_{n+1}{\parallel }^2+C\Delta t(\parallel {\sigma }_n{\parallel }^2+\parallel {\widehat{\sigma }}_n{\parallel }^2+\parallel {\widehat{\sigma }}_{n+1}{\parallel }^2).\hfill \end{array}$

Summing the inequality from 0 to $N-1$, we have

$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}\parallel {\sigma }_N{\parallel }^2+(1-C\epsilon )\Delta t\sum _{n=0}^{N-1}{\parallel \nabla {\sigma }_{n+{\displaystyle \frac{1}{2}}}\parallel }^2& \le C({(\Delta t)}^4+h^{2r+2})+C\parallel {\sigma }_0{\parallel }^2+C(\epsilon )\Delta t{\parallel {\sigma }_N\parallel }^2\hfill \\ \hfill & +C\Delta t\parallel {\widehat{\sigma }}_N{\parallel }^2+C\Delta t\sum _{n=0}^{N-1}(\parallel {\sigma }_n{\parallel }^2+\parallel {\widehat{\sigma }}_n{\parallel }^2).\hfill \end{array}$

Taking $\psi =\Delta t{\widehat{\sigma }}_{n+{\displaystyle \frac{1}{2}}}$ and using a similar discussion, we obtain

$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}\parallel {\widehat{\sigma }}_N{\parallel }^2+(1-C\epsilon )\Delta t\sum _{n=0}^{N-1}{\parallel \nabla {\widehat{\sigma }}_{n+{\displaystyle \frac{1}{2}}}\parallel }^2& \le C({(\Delta t)}^4+h^{2r+2})+C\parallel {\widehat{\sigma }}_0{\parallel }^2+C(\epsilon )\Delta t{\parallel {\widehat{\sigma }}_N\parallel }^2\hfill \\ \hfill & +C\Delta t\parallel {\sigma }_N{\parallel }^2+C\Delta t\sum _{n=0}^{N-1}(\parallel {\widehat{\sigma }}_n{\parallel }^2+\parallel {\sigma }_n{\parallel }^2).\hfill \end{array}$

Adding the two sides of the above inequalities, respectively, we can obtain

$\begin{array}{cc}\hfill & ({\displaystyle \frac{1}{2}}-C(\epsilon )\Delta t)(\parallel {\sigma }_N{\parallel }^2+\parallel {\widehat{\sigma }}_N{\parallel }^2)+(1-C\epsilon )\Delta t\sum _{n=0}^{N-1}(\parallel \nabla {\sigma }_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel \nabla {\widehat{\sigma }}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \\ \hfill & \le C({(\Delta t)}^4+h^{2r+2})+C(\parallel {\sigma }_0{\parallel }^2+\parallel {\widehat{\sigma }}_0{\parallel }^2)+C\Delta t\sum _{n=0}^{N-1}(\parallel {\widehat{\sigma }}_n{\parallel }^2+\parallel {\sigma }_n{\parallel }^2).\hfill \end{array}$

By taking ${\epsilon }_0,\Delta t$ to be sufficiently small such that $2\Delta t(1-C{\epsilon }_0)>0$$,{\displaystyle \frac{1}{2}}-C({\epsilon }_0)\Delta t>0$, and then using the discrete version of Gronwall’s inequality, we obtain

$\begin{array}{c}\hfill {∥{\sigma }_N∥}^2+{∥{\widehat{\sigma }}_N∥}^2\le C\left({\left(\Delta t\right)}^4+h^{2r+2}\right)+C\left({∥{\sigma }_0∥}^2+{∥{\widehat{\sigma }}_0∥}^2\right),\end{array}$

since ${∥{\sigma }_0∥}^2+{∥{\widehat{\sigma }}_0∥}^2\le C{h}^{2r+2}$, the above formula implies

(42)$\underset{0\le n\le N}{max}∥{\sigma }_n∥+∥{\widehat{\sigma }}_n∥\le C\left(h^{r+1}+{(\Delta t)}^2\right).$

Finally, using (42), (22), and the triangle inequality, we obtain the conclusion. □

5. L² and H¹ Estimates for Extrapolated Crank–Nicolson Galerkin Approximation

In this section, we present a modified version of the Crank–Nicolson Galerkin method from Section 4. The uniqueness of the approximation solution and the optimal-order error estimate, both in $L^2$ and $H^1$ norms, are obtained.

In the previous section, we obtain the $L^2$ optimal error estimates for Crank–Nicolson Galerkin approximations. To achieve more precise error estimates in the $H^1$ norm, we modify the numerical scheme accordingly. In the following discussion, we still use the notation in Section 4, and we define the generalized difference operator as $\mathcal{D}{\varphi }_n={\displaystyle \frac{3}{2}}{\varphi }_n-{\displaystyle \frac{1}{2}}{\varphi }_{n-1}$.

We obtain the extrapolated Crank–Nicolson Galerkin approximation of Problem (1) by replacing $U_{n+{\displaystyle \frac{1}{2}}},V_{n+{\displaystyle \frac{1}{2}}}$ with the extrapolated values from two previous time steps, $t_n$ and $t_{n-1}$, in the terms $g_1(u)u,g_2(v)v$, and the integral terms.

(43)$\left\{\begin{array}{cc}\left({\partial }_t{U}_n,\varphi \right)\hfill & +(\nabla U_{n+{\displaystyle \frac{1}{2}}},\nabla \varphi )+⟨g_1\left(\mathcal{D}{U}_n\right)\mathcal{D}{U}_n,\varphi ⟩\hfill \\ & =\left(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_1(\mathcal{D}{U}_n,\mathcal{D}{V}_n)\mathrm{d}x,\varphi \right),\forall \varphi \in S_h,\hfill \\ \left({\partial }_t{V}_n,\psi \right)\hfill & +(\nabla V_{n+{\displaystyle \frac{1}{2}}},\nabla \psi )+⟨g_2\left(\mathcal{D}{V}_n\right)\mathcal{D}{V}_n,\psi ⟩\hfill \\ & =\left(k_2(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_2(\mathcal{D}{U}_n,\mathcal{D}{V}_n)\mathrm{d}x,\psi \right),\forall \psi \in S_h.\hfill \end{array}\right.$

Since the above formula holds for $n\ge 1$, to obtain the optimal error estimates in both $L^2$ and $H^1$ norms, it is necessary to select $U_0,U_1,V_0,V_1$ to meet certain requirements.

We note that other extrapolation approaches, such as $2{\varphi }_{n-{\displaystyle \frac{1}{2}}}-{\varphi }_{n-{\displaystyle \frac{3}{2}}}$, exist. However, as shown by the Taylor expansions below, this alternative method results in a larger truncation error. Another drawback is that it involves three time node values rather than two, which increases the complexity.

$\begin{array}{cc}\hfill & {\displaystyle \frac{3}{2}}{\varphi }_n-{\displaystyle \frac{1}{2}}{\varphi }_{n-1}=\varphi (t_{n+{\displaystyle \frac{1}{2}}})-{\displaystyle \frac{3}{8}}{\varphi }_{tt}(t_{n+{\displaystyle \frac{1}{2}}}){(\Delta t)}^2+\mathcal{O}({(\Delta t)}^3),\hfill \\ \hfill & 2{\varphi }_{n-{\displaystyle \frac{1}{2}}}-{\varphi }_{n-{\displaystyle \frac{3}{2}}}=\varphi (t_{n+{\displaystyle \frac{1}{2}}})-{\displaystyle \frac{7}{8}}{\varphi }_{tt}(t_{n+{\displaystyle \frac{1}{2}}}){(\Delta t)}^2+\mathcal{O}({(\Delta t)}^3).\hfill \end{array}$

Theorem 7. 

The sequences ${\{U_n,V_n\}}_{n=0}^N$ defined by Equation (43) exist and are unique.

Proof. 

Firstly, we prove the existence of fully discrete solutions.

We use induction to prove the result. Assume that ${\left\{U_k,V_k\right\}}_{k=0}^n$ are given. For $\mathit{z}=\left(z_1,z_2\right)\in S_h\times S_h$, we define a mapping $F:\left(S_h\times S_h\right)\to \left(S_h\times S_h\right)$ such that

(44)$\begin{array}{cc}\hfill \left(\mathit{F}\left(\mathit{z}\right),\Psi \right)& =\left(\mathit{z}-{\mathit{x}}_{\mathit{n}},\Psi \right)+{\displaystyle \frac{1}{2}}\Delta t\left(\nabla \mathit{z},\nabla \Psi \right)+{\displaystyle \frac{1}{2}}\Delta t⟨\mathit{g}\left(\mathcal{D}{\mathit{x}}_{\mathit{n}}\right)\mathcal{D}{\mathit{x}}_{\mathit{n}},\Psi ⟩\hfill \\ & -{\displaystyle \frac{1}{2}}\Delta t\left(\mathit{k}{\int }_{\Omega }\mathit{f}\left(\mathcal{D}{U}_{\mathit{n}},\mathcal{D}{V}_{\mathit{n}}\right)\mathrm{d}x,\Psi \right)\hfill \end{array}$

for all $\Psi =\left(\varphi ,\psi \right)\in S_h\times S_h$, where ${\mathit{x}}_{\mathit{n}}=\left(U_n,V_n\right),\mathit{g}\left(\mathit{z}\right)\mathit{z}=\left(g_1\left(z_1\right)z_1,g_2\left(z_2\right)z_2\right)$, $\mathit{f}=\left(f_1\left(z_1,z_2\right),f_2\left(z_1,z_2\right)\right)$, and $\mathit{k}=\left(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}),k_2(x,t_{n+{\displaystyle \frac{1}{2}}})\right)$.

Taking $\Psi =\mathit{z}$ in the above formula, since $d=dim\Omega \le 3,H^1(\Omega )\to L^4(\partial \Omega )$ is a continuous embedding, we obtain

(45)$\begin{array}{cc}\hfill (F(\mathit{z}),\mathit{z})& =\left(\mathit{z}-{\mathit{x}}_n,\mathit{z}\right)+{\displaystyle \frac{1}{2}}\Delta t{\parallel \nabla \mathit{z}\parallel }^2+{\displaystyle \frac{1}{2}}\Delta t⟨\mathit{g}\left(\mathcal{D}{\mathit{x}}_n\right)\mathcal{D}{\mathit{x}}_n,\mathit{z}⟩-{\displaystyle \frac{1}{2}}\Delta t\left(\mathit{k}{\int }_{\Omega }\mathit{f}\left(\mathcal{D}{U}_n,\mathcal{D}{V}_n\right)\mathrm{d}x,\mathit{z}\right)\hfill \\ \hfill & \ge {\parallel \mathit{z}\parallel }^2-∥{\mathit{x}}_n∥\parallel \mathit{z}\parallel +{\displaystyle \frac{1}{2}}{\Delta t\parallel \nabla \mathit{z}\parallel }^2-{\displaystyle \frac{1}{4}}\Delta t{\left|\mathit{g}\left(\mathcal{D}{\mathit{x}}_n\right)\mathcal{D}{\mathit{x}}_n\right|}^2-{\displaystyle \frac{1}{4}}\Delta t{|\mathit{z}|}^2\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}\Delta t\parallel \mathit{k}(x,t_{n+{\displaystyle \frac{1}{2}}})\parallel ∥{\int }_{\Omega }\mathit{f}\left(z_1,z_2\right)\mathrm{d}x∥\parallel \mathit{z}\parallel \hfill \\ \hfill & \ge {\displaystyle \frac{1}{2}}\left({\parallel \mathit{z}\parallel }^2-{∥{\mathit{x}}_n∥}^2\right)+{\displaystyle \frac{1}{2}}\Delta t{\parallel \nabla \mathit{z}\parallel }^2-{\displaystyle \frac{1}{4}}\Delta t{\left|\mathit{g}\left(\mathcal{D}{\mathit{x}}_n\right)\mathcal{D}{\mathit{x}}_n\right|}^2-{\displaystyle \frac{1}{4}}\Delta t\left({\epsilon \parallel \nabla \mathit{z}\parallel }^2+C(\epsilon )\parallel \mathit{z}\parallel \right)\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}C\Delta t∥\mathit{f}\left(z_1,z_2\right)∥\parallel \mathit{z}\parallel \hfill \\ \hfill & \ge {\displaystyle \frac{1}{2}}\left({\parallel \mathit{z}\parallel }^2-{∥{\mathit{x}}_n∥}^2\right)+{\displaystyle \frac{1}{2}}\Delta t{\parallel \nabla \mathit{z}\parallel }^2-{\displaystyle \frac{1}{4}}\Delta t{\left|\mathit{g}\left(\mathcal{D}{\mathit{x}}_n\right)\mathcal{D}{\mathit{x}}_n\right|}^2-{\displaystyle \frac{1}{4}}\Delta t\left({\epsilon \parallel \nabla \mathit{z}\parallel }^2+C(\epsilon )\parallel \mathit{z}\parallel \right)\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}C\Delta t{\parallel \mathit{z}\parallel }^2-{\displaystyle \frac{1}{2}}C\Delta t\hfill \\ \hfill & \ge \left({\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{2}}C\Delta t-{\displaystyle \frac{1}{4}}C(\epsilon )\Delta t\right){\parallel \mathit{z}\parallel }^2+\left({\displaystyle \frac{1}{2}}\Delta t-{\displaystyle \frac{1}{4}}\epsilon \Delta t\right){\parallel \nabla \mathit{z}\parallel }^2-{\displaystyle \frac{1}{4}}\Delta t{\left|\mathit{g}\left(\mathcal{D}{\mathit{x}}_n\right)\mathcal{D}{\mathit{x}}_n\right|}^2\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{∥{\mathit{x}}_n∥}^2-{\displaystyle \frac{1}{2}}C\Delta t.\hfill \end{array}$

Taking ${\epsilon }_0,\Delta t$ to be sufficiently small such that $\left({\displaystyle \frac{1}{2}}\Delta t-{\displaystyle \frac{1}{4}}{\epsilon }_0\Delta t\right)>0,$$\left({\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{2}}C\Delta t-{\displaystyle \frac{1}{4}}C\left({\epsilon }_0\right)\Delta t\right)$$>{\displaystyle \frac{1}{4}},{\displaystyle \frac{1}{2}}C\Delta t<{\displaystyle \frac{1}{4}}$, and choosing ${∥\mathit{z}∥}^2>\Delta t{\left|\mathit{g}\left(\mathcal{D}{\mathit{x}}_{\mathit{n}}\right)\mathcal{D}{\mathit{x}}_{\mathit{n}}\right|}^2+2{∥{\mathit{x}}_{\mathit{n}}∥}^2+1$, we obtain

$\left(F\left(\mathit{z}\right),\mathit{z}\right)\ge {\displaystyle \frac{1}{4}}\left({∥\mathit{z}∥}^2-\Delta t{\left|\mathit{g}\left(\mathcal{D}{\mathit{x}}_{\mathit{n}}\right)\mathcal{D}{\mathit{x}}_{\mathit{n}}\right|}^2-2{∥{\mathit{x}}_{\mathit{n}}∥}^2-1\right)>0.$

By applying Brouwer’s fixed-point theorem, we derive that there exists ${\mathit{z}}^{*}=\left(z_1^{*},z_2^{*}\right)\in S_h\times S_h$ such that $F\left({\mathit{z}}^{*}\right)=0$. Let $U_{n+1}=2z_1^{*}-U_n,V_{n+1}=2z_2^{*}-V_n$, and it is easy to verify that $U_{n+1},V_{n+1}$ is a solution satisfying (29). Then, we generalize the assumptions to ensure the existence of the sequence of solutions.

We proceed with induction. First, we select $U_0={\widehat{U}}_0,U_1={\widehat{U}}_1,V_0={\widehat{V}}_0,V_1={\widehat{V}}_1$, such that $E_0={\widehat{E}}_0=E_1={\widehat{E}}_1=0$, and then we only need to prove $E_N={\widehat{E}}_N=0$ under the condition of $E_k={\widehat{E}}_k=0(\forall 0\le k<N)$.

Using the notation from Theorem 5, and following a similar method, we have

(46)$\begin{array}{cc}\hfill ({\partial }_t{E}_n,\varphi )& +(\nabla E_{n+{\displaystyle \frac{1}{2}}},\nabla \varphi )+⟨g_1(\mathcal{D}{U}_n)\mathcal{D}{U}_n-g_1(\mathcal{D}{\widehat{U}}_n)\mathcal{D}{\widehat{U}}_n,\varphi ⟩\hfill \\ \hfill & =\left(k_1(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_1(\mathcal{D}{U}_n,\mathcal{D}{V}_n)-f_1(\mathcal{D}{\widehat{U}}_n,\mathcal{D}{\widehat{V}}_n)\mathrm{d}x,\varphi \right),\hfill \end{array}$

(47)$\begin{array}{cc}\hfill ({\partial }_t{\widehat{E}}_n,\psi )& +(\nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}},\nabla \psi )+⟨g_2(\mathcal{D}{V}_n)\mathcal{D}{V}_n-g_2(\mathcal{D}{\widehat{V}}_n)\mathcal{D}{\widehat{V}}_n,\psi ⟩\hfill \\ \hfill & =\left(k_2(x,t_{n+{\displaystyle \frac{1}{2}}}){\int }_{\Omega }{f}_2(\mathcal{D}{U}_n,\mathcal{D}{V}_n)-f_2(\mathcal{D}{\widehat{U}}_n,\mathcal{D}{\widehat{V}}_n)\mathrm{d}x,\psi \right).\hfill \end{array}$

Let $\varphi =E_{k+{\displaystyle \frac{1}{2}}}$ in (46). Using a similar method as in Section 4, we have

(48)$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2\Delta t}}\left(\parallel E_{n+1}{\parallel }^2-{\parallel E_n\parallel }^2\right)+{\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill \le & C{\left|{\displaystyle \frac{3}{2}}{E}_n-{\displaystyle \frac{1}{2}}{E}_{n-1}\right|}^2+C{|E_{n+{\displaystyle \frac{1}{2}}}|}^2+C\left({∥{\displaystyle \frac{3}{2}}{E}_n-{\displaystyle \frac{1}{2}}{E}_{n-1}∥}^2+{∥{\displaystyle \frac{3}{2}}{\widehat{E}}_n-{\displaystyle \frac{1}{2}}{\widehat{E}}_{n-1}∥}^2+{\parallel E_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\right)\hfill \\ \hfill \le & C\epsilon (\parallel \nabla E_n{\parallel }^2+\parallel \nabla E_{n-1}{\parallel }^2+\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)+C(\epsilon )(\parallel E_n{\parallel }^2+\parallel E_{n-1}{\parallel }^2+\parallel E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \\ \hfill +& C(\parallel {\widehat{E}}_n{\parallel }^2+\parallel {\widehat{E}}_{n-1}{\parallel }^2)\hfill \\ \hfill \le & C(\parallel \nabla E_n{\parallel }^2+\parallel \nabla E_{n-1}{\parallel }^2+\parallel E_n{\parallel }^2+\parallel E_{n-1}{\parallel }^2+\parallel {\widehat{E}}_n{\parallel }^2+\parallel {\widehat{E}}_{n-1}{\parallel }^2)\hfill \\ \hfill +& C\epsilon \parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+C(\epsilon ){\parallel E_{n+1}\parallel }^2.\hfill \end{array}$

Multiplying (48) by $2\Delta t$ and summing from $n=1$ to $N-1$, we obtain

$\begin{array}{cc}\hfill (1-2C(\epsilon )\Delta t)\parallel E_N{\parallel }^2& +2\Delta t(1-C\epsilon )\sum _{n=1}^{N-1}{\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill & \le 2C\Delta t\sum _{n=1}^{N-1}\parallel \nabla E_n{\parallel }^2+2C\Delta t\sum _{n=1}^{N-1}(\parallel E_n{\parallel }^2+\parallel {\widehat{E}}_n{\parallel }^2).\hfill \end{array}$

Taking $\psi ={\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}$ in (47) and using a similar discussion, we obtain

$\begin{array}{cc}\hfill (1-2C(\epsilon )\Delta t)\parallel {\widehat{E}}_N{\parallel }^2& +2\Delta t(1-C\epsilon )\sum _{n=1}^{N-1}{\parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}\parallel }^2\hfill \\ \hfill & \le 2C\Delta t\sum _{n=1}^{N-1}\parallel \nabla {\widehat{E}}_n{\parallel }^2+2C\Delta t\sum _{n=1}^{N-1}(\parallel {\widehat{E}}_n{\parallel }^2+\parallel E_n{\parallel }^2).\hfill \end{array}$

Therefore, it is easy to obtain

$\begin{array}{cc}\hfill & (1-2C(\epsilon )\Delta t)(\parallel E_N{\parallel }^2+\parallel {\widehat{E}}_N{\parallel }^2)+2\Delta t(1-C\epsilon )\sum _{n=1}^{N-1}(\parallel \nabla E_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2+\parallel \nabla {\widehat{E}}_{n+{\displaystyle \frac{1}{2}}}{\parallel }^2)\hfill \\ \hfill & \le 2C\Delta t\sum _{n=1}^{N-1}(\parallel \nabla E_n{\parallel }^2+\parallel \nabla {\widehat{E}}_n{\parallel }^2)+2C\Delta t\sum _{n=1}^{N-1}(\parallel {\widehat{E}}_n{\parallel }^2+\parallel E_n{\parallel }^2).\hfill \end{array}$

Using a similar approach as in Section 4, and choosing sufficiently small $\epsilon ,\Delta t$ values, we ensure that the following inequality holds:

$\parallel E_N{\parallel }^2+{\parallel {\widehat{E}}_N\parallel }^2\le C\sum _{n=1}^{N-1}(\parallel \nabla E_n{\parallel }^2+\parallel \nabla {\widehat{E}}_n{\parallel }^2)+C\sum _{n=1}^{N-1}(\parallel {\widehat{E}}_n{\parallel }^2+\parallel E_n{\parallel }^2).$

Finally, the conclusion is established according to the inductive hypothesis. □

Theorem 8. 

Let ${\{U_n,V_n\}}_{n=0}^N$ be the approximate solution defined by (43). Assume that $U_0,U_1,V_0,V_1$ satisfy

(49) $\begin{array}{cc}\hfill \parallel U_0-{\mathcal{U}}_0{\parallel }_{H^1(\Omega )}+{\parallel U_1-{\mathcal{U}}_1\parallel }_{H^1(\Omega )}& \le C({(\Delta t)}^2+h^{r+1}),\hfill \\ \hfill \parallel V_0-{\mathcal{V}}_0{\parallel }_{H^1(\Omega )}+{\parallel V_1-{\mathcal{V}}_1\parallel }_{H^1(\Omega )}& \le C({(\Delta t)}^2+h^{r+1}).\hfill \end{array}$

Then, for a sufficiently small h and ${(\Delta t)}^4=o(h^d)$, the following error estimates hold:

$\begin{array}{cc}\hfill & \underset{2\le n\le N}{max}\parallel u_n-U_n\parallel +h\parallel u_n-U_n{\parallel }_{H^1(\Omega )}\le C(h^{r+1}+{(\Delta t)}^2),\hfill \\ \hfill & \underset{2\le n\le N}{max}\parallel v_n-V_n\parallel +h\parallel v_n-V_n{\parallel }_{H^1(\Omega )}\le C(h^{r+1}+{(\Delta t)}^2)\hfill \end{array}$

where C is a positive constant independent of h and $\Delta t$.

Proof. 

The conclusion of the theorem can be derived using a method similar to that in [12], with some differences in the details. Here, we outline the key steps rather than presenting the entire proof. First, from (1), (21) and (43), we can obtain

(50)$\begin{array}{cc}\hfill ({\partial }_t{\sigma }_k,\varphi )+(\nabla {\sigma }_{k+{\displaystyle \frac{1}{2}}},\nabla \varphi )& =-({\partial }_t{\rho }_k,\varphi )+\lambda ({\rho }_{k+{\displaystyle \frac{1}{2}}},\varphi )+({\partial }_t{u}_k-{(u_t)}_{k+{\displaystyle \frac{1}{2}}},\varphi )\hfill \\ & +⟨I_1^k,\varphi ⟩+⟨I_2^k,\varphi ⟩+⟨I_3^k,\varphi ⟩+\left(I_4^k,\varphi \right),\hfill \end{array}$

(51)$\begin{array}{cc}\hfill ({\partial }_t{\widehat{\sigma }}_k,\psi )+(\nabla {\widehat{\sigma }}_{k+{\displaystyle \frac{1}{2}}},\nabla \psi )& =-({\partial }_t{\widehat{\rho }}_k,\psi )+\lambda ({\widehat{\rho }}_{k+{\displaystyle \frac{1}{2}}},\psi )+({\partial }_t{v}_k-{(v_t)}_{k+{\displaystyle \frac{1}{2}}},\psi )\hfill \\ & +⟨J_1^k,\psi ⟩+⟨J_2^k,\psi ⟩+⟨J_3^k,\psi ⟩+\left(J_4^k,\varphi \right)\hfill \end{array}$

where

$\begin{array}{cc}\hfill I_1^k& =g_1\left(\mathcal{D}{U}_k\right)\mathcal{D}{U}_k-g_1\left(\mathcal{D}{\mathcal{U}}_k\right)\mathcal{D}{\mathcal{U}}_k,I_2^k=g_1\left(\mathcal{D}{\mathcal{U}}_k\right)\mathcal{D}{\mathcal{U}}_k-g_1\left({\mathcal{U}}_{k+{\displaystyle \frac{1}{2}}}\right){\mathcal{U}}_{k+{\displaystyle \frac{1}{2}}},\hfill \\ \hfill I_3^k& =g_1\left({\mathcal{U}}_{k+{\displaystyle \frac{1}{2}}}\right){\mathcal{U}}_{k+{\displaystyle \frac{1}{2}}}-{\left(g_1(\mathcal{U})\mathcal{U}\right)}_{k+{\displaystyle \frac{1}{2}}},J_1^k=g_2\left(\mathcal{D}{V}_k\right)\mathcal{D}{V}_k-g_2\left(\mathcal{D}{\mathcal{V}}_k\right)\mathcal{D}{\mathcal{V}}_k,\hfill \\ \hfill J_2^k& =g_2\left(\mathcal{D}{\mathcal{V}}_k\right)\mathcal{D}{\mathcal{V}}_k-g_2\left({\mathcal{V}}_{k+{\displaystyle \frac{1}{2}}}\right){\mathcal{V}}_{k+{\displaystyle \frac{1}{2}}},J_3^k=g_2\left({\mathcal{V}}_{k+{\displaystyle \frac{1}{2}}}\right){\mathcal{V}}_{k+{\displaystyle \frac{1}{2}}}-{\left(g_2(\mathcal{V})\mathcal{V}\right)}_{k+{\displaystyle \frac{1}{2}}},\hfill \\ \hfill I_4^k& ={\displaystyle \frac{1}{2}}\left(k_1\left(x,t_k\right){\int }_{\Omega }{f}_1\left(u_k,v_k\right)\mathrm{d}x+k_1\left(x,t_{k+1}\right){\int }_{\Omega }{f}_1\left(u_{k+1},v_{k+1}\right)\mathrm{d}x\right)\hfill \\ \hfill & -k_1\left(x,t_{k+{\displaystyle \frac{1}{2}}}\right){\int }_{\Omega }{f}_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)\mathrm{d}x,\hfill \end{array}$

$\begin{array}{cc}\hfill J_4^k& ={\displaystyle \frac{1}{2}}\left(k_2\left(x,t_k\right){\int }_{\Omega }{f}_2\left(u_k,v_k\right)\mathrm{d}x+k_2\left(x,t_{k+1}\right){\int }_{\Omega }{f}_2\left(u_{k+1},v_{k+1}\right)\mathrm{d}x\right)\hfill \\ \hfill & -k_2\left(x,t_{k+{\displaystyle \frac{1}{2}}}\right){\int }_{\Omega }{f}_2\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)\mathrm{d}x.\hfill \end{array}$

Note that

(52)$\begin{array}{cc}\hfill {∥{\partial }_t{u}_k-{\left({\displaystyle \frac{\partial u}{\partial t}}\right)}_{k+{\displaystyle \frac{1}{2}}}∥}^2& \le C{(\Delta t)}^3{\int }_{t_k}^{t_{k+1}}{∥u_{ttt}∥}^2\mathrm{d}s,\hfill \\ \hfill {∥{\partial }_t{v}_k-{\left({\displaystyle \frac{\partial v}{\partial t}}\right)}_{k+{\displaystyle \frac{1}{2}}}∥}^2& \le C{(\Delta t)}^3{\int }_{t_k}^{t_{k+1}}{∥v_{ttt}∥}^2\mathrm{d}s,\hfill \end{array}$

and

(53)$\Delta t{∥{\partial }_t{\rho }_k∥}^2\le {\int }_{t_k}^{t_{k+1}}{∥{\rho }_t∥}^2\mathrm{d}s,\Delta t{∥{\partial }_t{\widehat{\rho }}_k∥}^2\le {\int }_{t_k}^{t_{k+1}}{∥{\widehat{\rho }}_t∥}^2\mathrm{d}s.$

Taking $\varphi =D{\sigma }_k,\psi =D{\widehat{\sigma }}_k$ in (50) and (51), respectively, and summing from $k=1$ to n, we then use (52), (53), and (22) to obtain the following inequality:

$\begin{array}{cc}\hfill & {\displaystyle \frac{\Delta t}{2}}\sum _{k=1}^n{∥{\partial }_t{\sigma }_k∥}^2+{\displaystyle \frac{1}{2}}\sum _{k=1}^n\left({∥\nabla {\sigma }_{k+1}∥}^2-{∥\nabla {\sigma }_k∥}^2\right)\hfill \\ \hfill & \le C\left({(\Delta t)}^4+h^{2r+2}\right)+\left|\sum _{k=1}^n⟨I_1^k,D{\sigma }_k⟩\right|+\left|\sum _{k=1}^n⟨I_2^k,D{\sigma }_k⟩\right|+\left|\sum _{k=1}^n⟨I_3^k,D{\sigma }_k⟩\right|+\left|\sum _{k=1}^n\left(I_4^k,D{\sigma }_k\right)\right|\hfill \\ \hfill & =C\left({(\Delta t)}^4+h^{2r+2}\right)+\left|K_1\right|+\left|K_2\right|+\left|K_3\right|+\left|K_4\right|,\hfill \end{array}$

$\begin{array}{cc}\hfill & {\displaystyle \frac{\Delta t}{2}}\sum _{k=1}^n{∥{\partial }_t{\widehat{\sigma }}_k∥}^2+{\displaystyle \frac{1}{2}}\sum _{k=1}^n\left({∥\nabla {\widehat{\sigma }}_{k+1}∥}^2-{∥\nabla {\widehat{\sigma }}_k∥}^2\right)\hfill \\ \hfill & \le C\left({(\Delta t)}^4+h^{2r+2}\right)+\left|\sum _{k=1}^n⟨J_1^k,D{\widehat{\sigma }}_k⟩\right|+\left|\sum _{k=1}^n⟨J_2^k,D{\widehat{\sigma }}_k⟩\right|+\left|\sum _{k=1}^n⟨J_3^k,D{\widehat{\sigma }}_k⟩\right|+\left|\sum _{k=1}^n\left(J_4^k,D{\widehat{\sigma }}_k\right)\right|\hfill \\ \hfill & =C\left({(\Delta t)}^4+h^{2r+2}\right)+\left|L_1\right|+\left|L_2\right|+\left|L_3\right|+\left|L_4\right|.\hfill \end{array}$

Since the estimates of $K_1,K_2,K_3,L_1,L_2,L_3$ can be derived from [12], we only need to estimate $K_4,L_4$.

Note that

(54)$\begin{array}{cc}\hfill \left|K_4\right|& \le C\Delta t\sum _{k=1}^n{∥{\partial }_t{\sigma }_k∥}^2+C\Delta t\sum _{k=1}^n{∥M_1+M_2+M_3∥}^2\hfill \\ \hfill & \le C\Delta t\sum _{k=1}^n\left({∥{\partial }_t{\sigma }_k∥}^2+{∥M_1∥}^2+{∥M_2∥}^2+{∥M_3∥}^2\right),\hfill \end{array}$

where

$\begin{array}{cc}\hfill & M_1=\left[{\displaystyle \frac{k_1\left(x,t_k\right)+k_1\left(x,t_{k+1}\right)}{2}}-k_1\left(x,t_{k+{\displaystyle \frac{1}{2}}}\right)\right]{\int }_{\Omega }{f}_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)\mathrm{d}x,\hfill \\ \hfill & M_2=k_1\left(x,t_k\right){\int }_{\Omega }{\displaystyle \frac{f_1\left(u_k,v_k\right)+f_1\left(u_{k+1},v_{k+1}\right)}{2}}-f_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)\mathrm{d}x,\hfill \\ \hfill & M_3=\left(k_1\left(x,t_{k+1}\right)-k_1\left(x,t_k\right)\right){\int }_{\Omega }{f}_1\left(u_{k+1},v_{k+1}\right)-f_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)\mathrm{d}x.\hfill \end{array}$

Now, we estimate ${∥M_1∥}^2,{∥M_2∥}^2,{∥M_3∥}^2$, respectively. By using Taylor’s formula of $k_1(x,t_k),k_1(x,t_{k+1})$ at $t_{k+{\displaystyle \frac{1}{2}}}$, and eliminating the constant term and the linear term about $\Delta t$, we obtain

(55)$\begin{array}{cc}\hfill {∥M_1∥}^2& \le {∥{\int }_{\Omega }{f}_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)\mathrm{d}x∥}^2{∥{\displaystyle \frac{k_1\left(x,t_k\right)+k_1\left(x,t_{k+1}\right)}{2}}-k_1\left(x,t_{k+{\displaystyle \frac{1}{2}}}\right)∥}^2\hfill \\ \hfill & \le C{\left(\Delta t\right)}^4.\hfill \end{array}$

Similarly, we obtain

(56)$\begin{array}{cc}\hfill {∥M_2∥}^2& \le C{∥{\displaystyle \frac{f_1\left(u_k,v_k\right)+f_1\left(u_{k+1},v_{k+1}\right)}{2}}-f_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)∥}^2\hfill \\ \hfill & \le C{∥{\displaystyle \frac{f_1\left(u_k,v_k\right)+f_1\left(u_{k+1},v_{k+1}\right)}{2}}-f_1\left(u_{k+{\displaystyle \frac{1}{2}}},v_{k+{\displaystyle \frac{1}{2}}}\right)∥}^2\hfill \\ \hfill & +C{∥f_1\left(u_{k+{\displaystyle \frac{1}{2}}},v_{k+{\displaystyle \frac{1}{2}}}\right)-f_1\left(\mathcal{D}{u}_k,\mathcal{D}{v}_k\right)∥}^2+{∥f_1\left(\mathcal{D}{u}_k,\mathcal{D}{v}_k\right)-f_1\left(\mathcal{D}{\mathcal{U}}_k,\mathcal{D}{\mathcal{V}}_k\right)∥}^2\hfill \\ \hfill & +{∥f_1\left(\mathcal{D}{\mathcal{U}}_k,\mathcal{D}{\mathcal{V}}_k\right)-f_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)∥}^2\hfill \\ \hfill & \le C{(\Delta t)}^4+C{∥u_{k+{\displaystyle \frac{1}{2}}}-\mathcal{D}{u}_k∥}^2+C{∥v_{k+{\displaystyle \frac{1}{2}}}-\mathcal{D}{v}_k∥}^2\hfill \\ \hfill & +C{∥\mathcal{D}{\rho }_k∥}^2+C{∥\mathcal{D}{\widehat{\rho }}_k∥}^2+C{∥\mathcal{D}{\sigma }_k∥}^2+C{∥\mathcal{D}{\widehat{\sigma }}_k∥}^2\hfill \\ \hfill & \le C\left({(\Delta t)}^4+h^{2r+2}\right)+C{∥\mathcal{D}{\sigma }_k∥}^2+C{∥\mathcal{D}{\widehat{\sigma }}_k∥}^2,\hfill \end{array}$

and

(57)$\begin{array}{cc}\hfill {∥M_3∥}^2& \le C{∥\left(k_1\left(x,t_{k+1}\right)-k_1\left(x,t_k\right)\right)∥}^2{∥{\int }_{\Omega }{f}_1\left(u_{k+1},v_{k+1}\right)-f_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)\mathrm{d}x∥}^2\hfill \\ \hfill & \le C{(\Delta t)}^2{∥f_1\left(u_{k+1},v_{k+1}\right)-f_1\left(\mathcal{D}{U}_k,\mathcal{D}{V}_k\right)∥}^2\hfill \\ \hfill & \le C{(\Delta t)}^2\left({∥u_{k+1}-\mathcal{D}{u}_k∥}^2+{∥v_{k+1}-\mathcal{D}{v}_k∥}^2\right)\hfill \\ \hfill & +C{(\Delta t)}^2\left({∥\mathcal{D}{\rho }_k∥}^2+{∥\mathcal{D}{\widehat{\rho }}_k∥}^2+{∥\mathcal{D}{\sigma }_k∥}^2+{∥\mathcal{D}{\widehat{\sigma }}_k∥}^2\right)\hfill \\ \hfill & \le C\left({(\Delta t)}^4+h^{2r+2}\right)+C{(\Delta t)}^2\left({∥\mathcal{D}{\sigma }_k∥}^2+{∥\mathcal{D}{\widehat{\sigma }}_k∥}^2\right).\hfill \end{array}$

Using the Cauchy–Schwarz inequality and $\epsilon -$ inequality, we can easily obtain

(58)$\begin{array}{cc}\hfill \parallel {\sigma }_{n+1}{\parallel }^2& =\sum _{k=1}^n({\sigma }_{k+1}-{\sigma }_k,{\sigma }_{k+1}+{\sigma }_k)+{\parallel {\sigma }_1\parallel }^2\hfill \\ \hfill & \le {\left(\Delta t\sum _{k=1}^n{\parallel \partial {\sigma }_k\parallel }^2\right)}^{{\displaystyle \frac{1}{2}}}{\left(4\Delta t\sum _{k=1}^{n+1}{\parallel {\sigma }_k\parallel }^2\right)}^{{\displaystyle \frac{1}{2}}}+{\parallel {\sigma }_1\parallel }^2\hfill \\ \hfill & \le \epsilon \Delta t\sum _{k=1}^n\parallel \partial {\sigma }_k{\parallel }^2+C(\epsilon )\sum _{k=1}^{n+1}\Delta t\parallel {\sigma }_k{\parallel }^2+{\parallel {\sigma }_1\parallel }^2.\hfill \end{array}$

Then, by using (54)–(58), we have

$\begin{array}{cc}\hfill \left|K_4\right|& \le C({\left(\Delta t\right)}^4+h^{2r+2})+{\displaystyle \frac{\Delta t}{32}}\sum _{k=1}^n{∥{\partial }_t{\sigma }_k∥}^2+C\sum _{k=0}^n{s}_k{∥{\sigma }_k∥}^2+C\Delta t\sum _{k=1}^n\left({∥\mathcal{D}{\sigma }_k∥}^2+{∥\mathcal{D}{\widehat{\sigma }}_k∥}^2\right)\hfill \\ \hfill & \le C({\left(\Delta t\right)}^4+h^{2r+2})+{\displaystyle \frac{\Delta t}{16}}\sum _{k=1}^n{∥{\partial }_t{\sigma }_k∥}^2+C\sum _{k=0}^n{s}_k{∥{\sigma }_k∥}^2.\hfill \end{array}$

Here, $s_k$ is defined as

$s_k=\left\{\begin{array}{cc}& \Delta t,k=0,1,\hfill \\ \hfill & \Delta t+\sum _{j=k-2}^k\left({∥D{\sigma }_j∥}_{L^{\infty }(\Omega )}+{∥D{\rho }_j∥}_{L^{\infty }(\Omega )}+{∥D\nabla {\rho }_j∥}_{L^{\infty }(\Omega )}\right)\hfill \\ \hfill & +{\int }_{t_{k-2}}^{t_{k+1}}{∥{\mathcal{U}}_{ttt}∥}_{H^1(\Omega )}\mathrm{d}s,2\le k\le n-1,\hfill \\ \hfill & \Delta t+\sum _{j=n-2}^{n-1}{∥D{\sigma }_j∥}_{L^{\infty }(\Omega )}+\sum _{j=n-2}^n\left({∥D{\rho }_j∥}_{L^{\infty }(\Omega )}+{∥D\nabla {\rho }_j∥}_{L^{\infty }(\Omega )}\right)\hfill \\ \hfill & +{\int }_{t_{n-2}}^{t_{n+1}}{∥{\mathcal{U}}_{ttt}∥}_{H^1(\Omega )}\mathrm{d}s,k=n,\hfill \\ \hfill & \Delta t,k=n+1.\hfill \end{array}\right.$

Similarly, the estimation of $\left|L_4\right|$ can also be obtained. The remaining proof of the theorem is completely similar to the discussion in [12]. Therefore, it is omitted. □

6. Conclusions

Based on the discussion in this paper, the following conclusions can be drawn:

1. Under Conditions (A1)–(A3), the classical solution of System (1) has local existence and uniqueness in $C^{2+\alpha }\left({\overline{\Omega }}_T\right)\times C^{2+\alpha }\left({\overline{\Omega }}_T\right)$.

2. When the additional conditions $g_i\ge 0(i=1,2)$, $g_i\in C^3(\Omega )$ and (20) are satisfied, the three Galerkin approximation sequences of System (1) exist uniquely, and the following optimal-order error estimates hold:

   • (1). $\mathcal{O}(h^{r+1})$ in $L^2$ norm for continuous-time Galerkin approximation.

   • (2). $\mathcal{O}(h^{r+1}+{(\Delta t)}^2)$ in $L^2$ norm for the Crank–Nicolson Galerkin approximation.

   • (3). $\mathcal{O}(h^{r+1}+{(\Delta t)}^2)$ in both $L^2$ and $H^1$ norms for extrapolated Crank–Nicolson Galerkin approximation.

Author Contributions

Methodology, Q.G.; writing—original draft preparation, Y.Z.; formal analysis, B.Y. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

References

1. Gordon, P.V.; Ko, E.; Shivaji, R. Multiplicity and uniqueness of positive solutions for elliptic equations with nonlinear boundary conditions arising in a theory of thermal explosion. Nonlinear Anal. Real World Appl.; 2014; 15, pp. 51-57. [DOI: https://dx.doi.org/10.1016/j.nonrwa.2013.05.005]

2. Ko, E.; Prashanth, S. Positive solutions for elliptic equations in two dimensions arising in a theory of thermal explosion. Taiwan. J. Math.; 2015; 19, pp. 1759-1775. [DOI: https://dx.doi.org/10.11650/tjm.19.2015.5968]

3. Ko, E.; Ramaswamy, M.; Shivaji, R. Uniqueness of positive radial solutions for a class of semipositone problems on the exterior of a ball. J. Math. Anal. Appl.; 2015; 423, pp. 399-409. [DOI: https://dx.doi.org/10.1016/j.jmaa.2014.09.058]

4. Ma, W.; Yan, B. Global existence and uniform blow-up to a nonlocal parabolic system with nonlinear boundary conditions arising in a thermal explosion theory. Mathematics; 2023; 11, 1993. [DOI: https://dx.doi.org/10.3390/math11091993]

5. Ciarlet, P.G. The Finite Element Method for Elliptic Problems; SIAM: Philadelphia, PA, USA, 2002.

6. Baker, G.A. Error estimates for finite element methods for second order hyperbolic equations. SIAM J. Numer. Anal.; 1976; 13, pp. 564-576. [DOI: https://dx.doi.org/10.1137/0713048]

7. Li, Q.; Wei, H. The finite element method for nonlinear hyperbolic equation with nonlinear boundary condition. Math. Numer. Sin.; 1996; 18, pp. 285-294.

8. Shi, D.; Li, Z. Superconvergence analysis of the finite element method for nonlinear hyperbolic equations with nonlinear boundary condition. Appl. Math.-A J. Chin. Univ.; 2008; 23, pp. 455-462. [DOI: https://dx.doi.org/10.1007/s11766-008-1901-6]

9. Finlayson, B.A. Convergence of the Galerkin method for nonlinear problems involving chemical reaction. SIAM J. Numer. Anal.; 1971; 8, pp. 316-324. [DOI: https://dx.doi.org/10.1137/0708032]

10. Rachford, H., Jr. Two-level discrete-time Galerkin approximations for second order nonlinear parabolic partial differential equations. SIAM J. Numer. Anal.; 1973; 10, pp. 1010-1026. [DOI: https://dx.doi.org/10.1137/0710084]

11. Douglas, J.; Dupont, T. Galerkin methods for parabolic equations with nonlinear boundary conditions. Numer. Math.; 1973; 20, pp. 213-237. [DOI: https://dx.doi.org/10.1007/BF01436565]

12. Luskin, M. A Galerkin method for nonlinear parabolic equations with nonlinear boundary conditions. SIAM J. Numer. Anal.; 1979; 16, pp. 284-299. [DOI: https://dx.doi.org/10.1137/0716021]

13. Lin, Y. Galerkin methods for nonlinear parabolic integrodifferential equations with nonlinear boundary conditions. SIAM J. Numer. Anal.; 1990; 27, pp. 608-621. [DOI: https://dx.doi.org/10.1137/0727037]

14. Bu, W.; Tang, Y.; Yang, J. Galerkin finite element method for two-dimensional Riesz space fractional diffusion equations. J. Comput. Phys.; 2014; 276, pp. 26-38. [DOI: https://dx.doi.org/10.1016/j.jcp.2014.07.023]

15. Nedaiasl, K.; Dehbozorgi, R. Galerkin finite element method for nonlinear fractional differential equations. Numer. Algorithms; 2021; 88, pp. 113-141. [DOI: https://dx.doi.org/10.1007/s11075-020-01032-2]

16. Babuska, I.; Tempone, R.; Zouraris, G.E. Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM J. Numer. Anal.; 2004; 42, pp. 800-825. [DOI: https://dx.doi.org/10.1137/S0036142902418680]

17. Yan, Y. Galerkin finite element methods for stochastic parabolic partial differential equations. SIAM J. Numer. Anal.; 2005; 43, pp. 1363-1384. [DOI: https://dx.doi.org/10.1137/040605278]

18. Cheng, Y.; Shu, C.W. A discontinuous Galerkin finite element method for time dependent partial differential equations with higher order derivatives. Math. Comput.; 2008; 77, pp. 699-730. [DOI: https://dx.doi.org/10.1090/S0025-5718-07-02045-5]

19. Li, C.; Li, Z.; Wang, Z. Mathematical analysis and the local discontinuous Galerkin method for Caputo–Hadamard fractional partial differential equation. J. Sci. Comput.; 2020; 85, pp. 1-27. [DOI: https://dx.doi.org/10.1007/s10915-019-01102-1]

20. Pani, A.K. An H1-Galerkin mixed finite element method for parabolic partial differential equations. SIAM J. Numer. Anal.; 1998; 35, pp. 712-727. [DOI: https://dx.doi.org/10.1137/S0036142995280808]

21. Zhou, Z. An H1-Galerkin mixed finite element method for a class of heat transport equations. Appl. Math. Model.; 2010; 34, pp. 2414-2425. [DOI: https://dx.doi.org/10.1016/j.apm.2009.11.007]

22. Abbaszadeh, M.; Dehghan, M.; Zhou, Y. Crank–Nicolson/Galerkin spectral method for solving two-dimensional time-space distributed-order weakly singular integro-partial differential equation. J. Comput. Appl. Math.; 2020; 374, 112739. [DOI: https://dx.doi.org/10.1016/j.cam.2020.112739]

23. Xie, Z.; Li, X.; Tang, T. Convergence analysis of spectral Galerkin methods for Volterra type integral equations. J. Sci. Comput.; 2012; 53, pp. 414-434. [DOI: https://dx.doi.org/10.1007/s10915-012-9577-8]

24. Djomegne, L.; Kenne, C. Stackelberg exact controllability of a class of nonlocal parabolic equations. ESAIM Control. Optim. Calc. Var.; 2024; 30, 57. [DOI: https://dx.doi.org/10.1051/cocv/2024045]

25. Chen, L.Y.; Wu, H.Y.; Jiang, L.H. Ring-like two-breather structures of a partially nonlocal NLS system with different two-directional diffractions under a parabolic potential. Chaos Solitons Fractals; 2024; 178, 114330. [DOI: https://dx.doi.org/10.1016/j.chaos.2023.114330]

26. Liu, M.; Huang, P.; He, Y. A linearized Crank–Nicolson/Leapfrog scheme for the Landau–Lifshitz equation. Rocky Mt. J. Math.; 2023; 53, pp. 821-837. [DOI: https://dx.doi.org/10.1216/rmj.2023.53.821]

27. Wang, P.; Huang, P. Convergence of the Crank-Nicolson extrapolation scheme for the Korteweg-de Vries equation. Appl. Numer. Math.; 2019; 143, pp. 88-96. [DOI: https://dx.doi.org/10.1016/j.apnum.2019.04.001]

28. Atmani, S.; Biroud, K.; Daoud, M.; Laamri, E.H. On some fractional parabolic reaction-diffusion systems with gradient source terms. Fract. Calc. Appl. Anal.; 2024; 27, pp. 2644-2687. [DOI: https://dx.doi.org/10.1007/s13540-024-00316-x]

29. Hu, X.; Huang, P.; Feng, X. A new mixed finite element method based on the Crank-Nicolson scheme for Burgers’ equation. Appl. Math.; 2016; 61, pp. 27-45. [DOI: https://dx.doi.org/10.1007/s10492-016-0120-3]

30. Evans, L.C. Partial Differential Equations; American Mathematical Society: Providence, RI, USA, 2022; Volume 19.

31. Gómez, J.L.; Márquez, V.; Wolanski, N. Blow up results and localization of blow up points for the heat equation with a nonlinear boundary condition. J. Differ. Equ.; 1991; 92, pp. 384-401. [DOI: https://dx.doi.org/10.1016/0022-0396(91)90056-F]

32. Lieberman, G.M. Second Order Parabolic Differential Equations; World Scientific: Singapore, 1996.

33. Gilbarg, D.; Trudinger, N.S.; Gilbarg, D.; Trudinger, N. Elliptic Partial Differential Equations of Second Order; Springer: Berlin/Heidelberg, Germany, 1977; Volume 224.