# Finite Element Analysis of Time-Fractional Integro-differential Equation of
Kirchhoff type for Non-homogeneous Materials

Lalit Kumar 111Department of Mathematics, Indian Institute of Technology
Bombay, Mumbai-400076, India<EMAIL_ADDRESS>Sivaji Ganesh
Sista222Department of Mathematics, Indian Institute of Technology Bombay,
Mumbai-400076, India<EMAIL_ADDRESS>Konijeti Sreenadh††footnotemark:
333Department of Mathematics, Indian Institute of Technology Delhi, New
Delhi-110016, India<EMAIL_ADDRESS>
###### Abstract

In this paper, we study a time-fractional initial-boundary value problem of
Kirchhoff type involving memory term for non-homogeneous materials
($D_{\alpha}$). The energy argument is applied to derive the a priori bounds
on the solution of the problem ($D_{\alpha}$). Consequently, we prove the
existence and uniqueness of the weak solution to the considered problem. We
keep the time variable continuous and discretize the space domain using a
conforming FEM to obtain the semi discrete formulation of the problem under
consideration. The semi discrete error analysis is carried out by modifying
the standard Ritz-Volterra projection operator. To obtain the numerical
solution of the problem ($D_{\alpha}$) efficiently, we develop a new
linearized L1 Galerkin FEM. This numerical scheme is shown to have a
convergence rate of $O(h+k^{2-\alpha})$, where $\alpha~{}(0<\alpha<1)$ is the
fractional derivative exponent, $h$ and $k$ are the discretization parameters
in the space and time directions respectively. Further, this convergence rate
is improved in the time direction by proposing a novel linearized L2-1σ
Galerkin FEM. We prove that this numerical scheme has an accuracy rate of
$O(h+k^{2})$. Finally, a numerical experiment is conducted to validate our
theoretical claims.

Keywords: Nonlocal, Finite element method (FEM), Fractional time derivative,
Fractional Crank-Nicolson scheme, Integro-differential equation.
AMS subject classification. 34K30, 26A33, 65R10, 60K50.

##  1 Introduction

Let $\Omega$ be a convex and bounded subset of $\mathbb{R}^{d}~{}(d\geq 1)$
with smooth boundary $\partial\Omega$ and $[0,T]$ is a fixed finite time
interval. We consider the following integro-differential equation of Kirchhoff
type involving fractional time derivative of order $\alpha~{}(0<\alpha<1)$ for
non-homogeneous materials

$\displaystyle^{C}D^{\alpha}_{t}u-\nabla\cdot\left(M\left(x,t,\|\nabla
u\|^{2}\right)\nabla
u\right)=f(x,t)+\int_{0}^{t}b(x,t,s)u(s)~{}ds\quad\text{in}~{}\Omega\times(0,T],$
($D_{\alpha}$)

with initial and boundary conditions

$\displaystyle u(x,0)$ $\displaystyle=u_{0}(x)\quad\text{in}~{}\Omega,$
$\displaystyle u(x,t)$
$\displaystyle=0\quad\text{on}~{}\partial\Omega\times[0,T],$

where $u:=u(x,t):\Omega\times[0,T]\rightarrow\mathbb{R}$ is the unknown
function, $M:\bar{\Omega}\times[0,T]\times(0,\infty)\rightarrow(0,\infty),$
initial data $u_{0}$, source term $f$ are known functions and $b(x,t,s)$ is a
memory operator to be defined in Section 2. The notation
${}^{C}D_{t}^{\alpha}u$ in the problem ($D_{\alpha}$) is the fractional time
derivative of order $\alpha$ in the Caputo sense, which is defined in [30] as

$^{C}D_{t}^{\alpha}u=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{1}{(t-s)^{\alpha}}\frac{\partial
u}{\partial s}(s)~{}ds,$ (1.1)

where $\Gamma(\cdot)$ denotes the gamma function.

For the case $\alpha=1$, authors in [16] proposed a linearized backward Euler-
Galerkin FEM and a linearized Crank-Nicolson-Galerkin FEM with an accuracy
rate of $O(h+k)$ and $O(h+k^{2})$, respectively.

There are various notions of fractional derivatives other than the Caputo
derivative, which include Riemann-Liouville, Grünwald-Letnikov, Weyl,
Marchaud, and Riesz fractional derivatives [26, 30]. Among these, Caputo
fractional derivative and Riemann-Liouville fractional derivative are the most
commonly used in the literature [30]. These two fractional derivatives are
related to each other by the following relation for absolutely continuous
function $u$, see [30]

$^{C}D_{t}^{\alpha}u(t)=~{}^{R}D_{t}^{\alpha}\left(u(t)-u(0)\right):=\partial^{\alpha}_{t}u(t),$
(1.2)

where ${}^{R}D_{t}^{\alpha}u$ is the Riemann-Liouville fractional derivative
defined by

${}^{R}D_{t}^{\alpha}u(t)=\frac{1}{\Gamma(1-\alpha)}\frac{d}{dt}\int_{0}^{t}\frac{1}{(t-s)^{\alpha}}u(s)~{}ds.$
(1.3)

It is worth noting that fractional derivatives other than the Caputo
fractional derivative require initial condition containing the limiting value
of fractional derivative at $t=0$ [30], which has no physical interpretation.
An advantage of choosing Caputo fractional derivative in the problem
($D_{\alpha}$) is that it allows the initial and boundary conditions in the
same way as those for integer order differential equations.

There are many physical and biological processes in which the mean-squared
displacement of the particle motion grows only sublinearly with time $t$,
instead of linear growth. For instance, acoustic wave propagation in
viscoelastic materials [22], cancer invasion system [23], anomalous diffusion
transport [25], which cannot be described accurately by classical models
having integer order derivatives. Therefore the study of fractional
differential equations has evolved immensely in recent years.

Mathematical problems involving fractional time derivatives have been studied
by many researchers, for instance, see [9, 13, 31]. Analytical solutions of
fractional differential equations are expressed in terms of Mittag-Leffler
function, Fox $H$-functions, Green functions, and hypergeometric functions.
Such special functions are more complex to compute, which restrict the
applications of fractional calculus in applied sciences. This motivates the
researchers to develop numerical algorithms for solving fractional
differential equations.

There are two predominant discretization techniques in time for fractional
differential equations, namely Grünwald-Letnikov approximation [8] and L1 type
approximation [1, 20]. The second category, viz. the L1 type approximation
scheme is based on piecewise interpolation of the integrand in the definition
(1.1) of the Caputo fractional derivative. Lin and Xu in [20] developed the L1
scheme based on piecewise linear interpolation for Caputo fractional
derivative and Legendre spectral method in space for the following time-
fractional PDE in one space dimension

$\displaystyle{}^{C}D^{\alpha}_{t}u-\frac{\partial^{2}u}{\partial x^{2}}$
$\displaystyle=f(x,t)\quad x\in(0,1),~{}t\in(0,T],$ (1.4) $\displaystyle
u(x,0)$ $\displaystyle=g(x)\quad x\in(0,1),$ $\displaystyle u(0,t)=u(1,t)$
$\displaystyle=0\quad 0\leq t\leq T,$

and achieved the convergence estimates of $O(h^{2}+k^{2-\alpha})$ for
solutions belonging to $C^{2}\left([0,T];H^{2}(\Omega)\cap
H^{1}_{0}(\Omega)\right)$. Recently, Alikhanov [1] proposed a modification of
the L1 type scheme in the time direction and difference scheme in the space
direction for some linear extension the problem (1.4). In his work, the author
proved that the convergence rate is of $O(h^{2}+k^{2})$ for solutions
belonging to $C^{3}\left([0,T];C^{4}(\Omega)\right)$.

On a similar note, there has been considerable attention devoted to the
nonlocal diffusion problems where diffusion coefficient depends on the entire
domain rather than pointwise. Lions [21] studied the following problem

$\frac{\partial^{2}u}{\partial t^{2}}-M\left(x,\int_{\Omega}|\nabla
u|^{2}dx\right)\Delta u=f(x,t)\quad\text{in}\quad\Omega\times[0,T],$

which models transversal oscillations of an elastic string or membrane by
considering the change in length during vibrations. Also, the nonlocal terms
appear in various physical and biological systems. For instance, the
temperature in a thin region during friction welding [14], Kirchhoff equations
with magnetic field [27], Ohmic heating with variable thermal conductivity
[33], and many more. We cite [6, 10, 28] for some contemporary works related
to the existence, uniqueness, and regularity of the problems involving
Kirchhoff type diffusion coefficient.

The models discussed above behave accurately only for a perfectly homogeneous
medium, but in real-life situations, a large number of heterogeneities are
present, which cause some memory effect or feedback term. These phenomena
cannot be described by classical PDEs, which motivate us to study time-
fractional PDEs for non-homogeneous materials. We note that this class of
equations has not been analyzed in the literature yet, and this is the first
attempt to establish new results for the problem ($D_{\alpha}$).

The salient features of the considered problem is its doubly nonlocal nature
due to the presence of Kirchhoff term and fractional derivative or memory
term. The appearance of the Kirchhoff term creates major difficulties in the
theoretical as well as in numerical analysis of the problem under
consideration. We cannot apply Laplace/Fourier transformation in the problem
($D_{\alpha}$), therefore explicit representation of its solution in terms of
Fourier expansion is not possible. To resemble this issue, we use Galerkin
method to show the well-posedness of the weak formulation of the problem
($D_{\alpha}$).

The fully discrete formulation of the considered problem produces a system of
nonlinear algebraic equations. In general, numerical schemes based on the
Newton method are adopted to solve this system [12]. The Kirchhoff term leads
to the highly non-sparse Jacobian of this system [12]. As a result of which we
require high computational cost as well as huge computer storage for solving
this system. We reduce these costs by developing a new linearization technique
for the nonlinearity.

On the other hand, the memory term incorporates the history of the phenomena
under investigation by virtue of which we need to store the value of
approximate solution at all previous time steps. This process demands large
computer memory. We overcome this difficulty by discretizing the memory term
using modified trapezoidal or modified Simpson’s rule [29].

To prove the well-posedness of the weak formulation of the problem
($D_{\alpha}$), we reduce the weak formulation onto a finite dimensional
subspace of $H^{1}_{0}(\Omega)$. The theory of fractional differential
equations [7] ensures the existence of Galerkin sequence of weak solutions.
The a priori bounds on these Galerkin sequences are attained by employing the
energy argument. We make use of these a priori bounds in Aubin-Lions type
compactness lemma [17] to prove that the Galerkin sequence converges to the
weak solution of the problem ($D_{\alpha}$).

To determine the semi discrete error estimates, we alter the definition of
standard Ritz-Volterra projection operator [4] so that it reduces the
complications caused by the Kirchhoff term. This modified Ritz-Volterra
projection operator follows the best approximation properties same as that of
standard one. These best approximation properties play a key role in deriving
the semi discrete error estimates.

To obtain the numerical solution, we construct two fully discrete formulations
for the problem ($D_{\alpha}$) by discretizing the space domain using a
conforming FEM [32] and the time direction by uniform mesh. First, we develop
a new linearized L1 Galerkin FEM. This method comprises of L1 type
approximation [20] of the Caputo fractional derivative, linearization
technique for the Kirchhoff type nonlinearity, and modified Simpson’s rule
[29] for approximation of the memory term. We acquire the a priori bounds on
the solution of this numerical scheme and show that this numerical scheme is
accurate of $O(h+k^{2-\alpha})$.

Further, we increase the accuracy of this scheme in the time direction by
replacing the L1 scheme with the L2-1σ scheme [1] for the approximation of the
Caputo fractional derivative. As a consequence, we propose a new linearized
L2-1σ Galerkin FEM which has a convergence rate of $O(h+k^{2})$. These
numerical results are supported by conducting a numerical experiment in MATLAB
software.

Turning to the layout of this paper: In Section 2, we provide some notations,
assumptions, and preliminaries results that will be used throughout this work.
In Section 3, we state main contributions of this article. Section 4 contains
the proof of well-posedness of the weak formulation of the problem
($D_{\alpha}$). In Section 5, we define semi discrete formulation of the
considered problem and derive a priori bounds as well as error estimates on
semi discrete solutions. In Section 6, we develop a new linearized L1 Galerkin
FEM. We derive a priori bounds on numerical solutions of the developed
numerical scheme and prove its accuracy rate of $O(h+k^{2-\alpha})$. In
Section 7, we achieve improved convergence rate of $O(h+k^{2})$ by proposing a
new linearized L2-1σ Galerkin FEM. Section 8 includes a numerical experiment
that confirms the sharpness of theoretical results. Finally, we conclude this
work in Section 9.

##  2 Preliminaries

Let $L^{1}(\Omega)$ be the set of all equivalence classes of the integrable
functions on $\Omega$ with the norm

$\|g\|_{L^{1}(\Omega)}=\int_{\Omega}|g(x)|~{}dx~{}~{}\text{for}~{}~{}g\in
L^{1}(\Omega).$ (2.1)

Let $L^{2}(\Omega)$ be the set of all equivalence classes of the square
integrable functions on $\Omega$ with the norm

$\|g\|^{2}=\int_{\Omega}|g(x)|^{2}~{}dx~{}~{}\text{for}~{}~{}g\in
L^{2}(\Omega).$ (2.2)

The norm defined in (2.2) is induced by the inner product $(\cdot,\cdot)$ as
follows

$(g,h)=\int_{\Omega}g(x)h(x)~{}dx~{}~{}\text{for}~{}~{}g,h\in L^{2}(\Omega).$
(2.3)

The Sobolev space $W^{1,1}(\Omega)$ is the collection of all functions in
$L^{1}(\Omega)$ such that its distributional derivative of order one is also
in $L^{1}(\Omega)$, i.e.,

$W^{1,1}(\Omega)=\left\\{g\in L^{1}(\Omega);~{}Dg\in L^{1}(\Omega)\right\\}.$
(2.4)

The norm on the space $W^{1,1}(\Omega)$ is given by

$\|g\|_{W^{1,1}(\Omega)}=\|g\|_{L^{1}(\Omega)}+\|Dg\|_{L^{1}(\Omega)}~{}~{}\text{for}~{}~{}g\in
W^{1,1}(\Omega).$ (2.5)

The sobolev space $H^{m}(\Omega),\left(m\in\\{1,2\\}\right)$ is the set of all
functions in $L^{2}(\Omega)$ such that its distributional derivatives upto
order $m$ are also in $L^{2}(\Omega)$, i.e.,

$H^{m}(\Omega)=\left\\{g\in L^{2}(\Omega);~{}D^{\beta}g\in
L^{2}(\Omega),~{}|\beta|\leq m\right\\},$ (2.6)

where $\beta$ is multiindex. The norm on the space $H^{m}(\Omega)$ is induced
by the following inner product $(\cdot,\cdot)_{m}$ as follows

$(g,h)_{m}=\sum_{|\beta|\leq
m}(D^{\beta}g,D^{\beta}h)~{}~{}\text{for}~{}~{}g,h\in H^{m}(\Omega).$ (2.7)

We denote $H^{m}_{0}(\Omega),\left(m\in\\{1,2\\}\right)$ be the closure of
$C^{\infty}_{C}(\Omega)$ in $H^{m}(\Omega)$. The space $H^{m}_{0}(\Omega)$ can
be characterised by the functions in $H^{m}(\Omega)$ having zero trace on the
boundary $\partial\Omega$ [15, Section 2.7]. The dual space of the
$H^{m}_{0}(\Omega)$ is denoted by $H^{-m}(\Omega)$.

For any Hilbert space $X$, we denote $L^{2}(0,T;X)$ be the set of all
measurable functions $g:[0,T]\rightarrow X$ such that

$\int_{0}^{T}\|g(s)\|^{2}_{X}~{}ds<\infty.$ (2.8)

The norm on the space $L^{2}(0,T;X)$ is given by

$\|g\|^{2}_{L^{2}(0,T;X)}=\int_{0}^{T}\|g(s)\|_{X}^{2}~{}ds~{}~{}\text{for}~{}~{}g\in
L^{2}(0,T;X).$ (2.9)

We also define a weighted $L^{2}_{\alpha}(0,T;X)$ space consisting of all
measurable functions $g:[0,T]\rightarrow X$ such that

$\sup_{t\in(0,T)}\left(\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}\|g(s)\|^{2}_{X}~{}ds\right)<\infty.$
(2.10)

The norm on the space $L^{2}_{\alpha}(0,T;X)$ is given by [18, (4.5)]

$\|g\|^{2}_{L^{2}_{\alpha}(0,T;X)}=\sup_{t\in(0,T)}\left(\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}\|g(s)\|^{2}_{X}~{}ds\right)~{}~{}\text{for}~{}~{}g\in
L^{2}_{\alpha}(0,T;X).$ (2.11)

One can observe that $L^{2}_{\alpha}(0,T;X)\subset L^{2}(0,T;X)$. The set of
all measurable functions $g:[0,T]\rightarrow X$ such that

$\text{ess}\sup_{t\in(0,T)}\|g(t)\|_{X}<\infty$ (2.12)

is denoted by $L^{\infty}(0,T;X)$. The norm on this space is given by

$\|g\|_{L^{\infty}(0,T;X)}=\text{ess}\sup_{t\in(0,T)}\|g(t)\|_{X}~{}~{}\text{for}~{}~{}g\in
L^{\infty}(0,T;X).$ (2.13)

For any two quantities $a$ and $b$, the notation $a\lesssim b$ means that
there exists a generic positive constant $C$ such that $a\leq Cb$, where $C$
depends on data but independent of discretization parameters and may vary at
different occurrences.

Throughout the paper, we assume the following hypotheses on data:
(H1) Initial data $u_{0}\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ and source
term $f\in L^{\infty}\left(0,T;L^{2}(\Omega)\right)$.
(H2) Diffusion coefficient
$M:\bar{\Omega}\times[0,T]\times(0,\infty)\rightarrow(0,\infty)$ is a
Lipschitz continuous function such that there exists a constant $m_{0}$ which
satisfies

$M(x,t,s)\geq m_{0}>0~{}\text{for all
}~{}(x,t,s)\in\bar{\Omega}\times[0,T]\times(0,\infty)~{}\text{and}\left(m_{0}-4L_{M}K^{2}\right)>0,$

where $K=\left(\|\nabla
u_{0}\|+\|f\|_{L^{\infty}\left(0,T;L^{2}(\Omega)\right)}\right)$ and $L_{M}$
is a Lipschitz constant.
(H3) Memory operator $b(x,t,s)$ is a second order partial differential
operator of the form

$b(x,t,s)u(s):=-\nabla\cdot(b_{2}(x,t,s)\nabla
u(s))+\nabla\cdot(b_{1}(x,t,s)u(s))+b_{0}(x,t,s)u(s),$

with $b_{2}:\bar{\Omega}\times[0,T]\times[0,T]\rightarrow\mathbb{R}^{d\times
d}$ is a symmetric and positive definite matrix with entries
$[b_{2}^{ij}(x,t,s)]$,
$b_{1}:\bar{\Omega}\times[0,T]\times[0,T]\rightarrow\mathbb{R}^{d}$ is a
vector with entries $[b_{1}^{j}(x,t,s)]$ and
$b_{0}:\bar{\Omega}\times[0,T]\times[0,T]\rightarrow\mathbb{R}$ is a scalar
function. We assume that $b_{2}^{ij},b_{1}^{j},b_{0}$ are smooth functions in
all variables $(x,t,s)\in\bar{\Omega}\times[0,T]\times[0,T]$ for
$i,j=1,2,\dots,d$.

We define a function $B(t,s,u(s),v)$ for all $t,s$ in $[0,T]$ and for all
$u(s),v$ in $H^{1}_{0}(\Omega)$ as follows

$B(t,s,u(s),v):=(b_{2}(x,t,s)\nabla u(s),\nabla
v)+(\nabla\cdot(b_{1}(x,t,s)u(s)),v)+(b_{0}(x,t,s)u(s),v).$ (2.14)

Using (H3) and Poincaré inequality we can prove that there exists a positive
constant $B_{0}$ such that for all $t,s$ in $[0,T]$ and for all $u,v$ in
$H^{1}_{0}(\Omega)$ we have

$|B(t,s,u(s),v)|\leq B_{0}\|\nabla u(s)\|~{}\|\nabla v\|.$ (2.15)

We denote

$k(t):=\frac{t^{-\alpha}}{\Gamma(1-\alpha)},$ (2.16)

and $\ast$ indicates the convolution of two integrable functions $g$ and $h$
on $[0,T]$ as

$(g\ast h)(t)=\int_{0}^{t}g(t-s)h(s)~{}ds\quad\text{for
all}~{}~{}t~{}~{}\text{in}~{}~{}[0,T].$ (2.17)

######  Remark 1.

Note that $l(t)$ defined by $l(t):=\frac{t^{\alpha-1}}{\Gamma(\alpha)}$
satisfies $k\ast l=1$.

######  Lemma 2.1.

[11, Lemma 18.4.1] Let $H$ be a real Hilbert space and $T>0$. Then for any
$\tilde{k}\in W^{1,1}(0,T)$ and $v\in L^{2}(0,T;H)$ we have

$\displaystyle\left(\frac{d}{dt}\left(\tilde{k}\ast
v\right)(t),v(t)\right)_{H}$
$\displaystyle=\frac{1}{2}\frac{d}{dt}\left(\tilde{k}\ast\|v\|^{2}_{H}\right)(t)+\frac{1}{2}\tilde{k}(t)\|v\|_{H}^{2}$
(2.18)
$\displaystyle+\frac{1}{2}\int_{0}^{t}\left[-\tilde{k}^{\prime}(s)\right]\|v(t)-v(t-s)\|_{H}^{2}~{}ds~{}~{}\text{a.e.}~{}~{}t\in(0,T).$

######  Lemma 2.2.

[2, Theorem 8] Let $u,v$ be two nonnegative integrable functions on $[a,b]$
and $g$ a continuous function in $[a,b]$. Assume that $v$ is nondecreasing in
$[a,b]$ and $g$ is nonnegative and nondecreasing in $[a,b]$. If

$u(t)\leq
v(t)+g(t)\int_{a}^{t}(t-s)^{\alpha-1}u(s)~{}ds~{}~{}\text{for}~{}\alpha\in(0,1)~{}\text{and}~{}\forall~{}t\in[a,b],$

then

$u(t)\leq
v(t)E_{\alpha}\left[g(t)\Gamma(\alpha)(t-a)^{\alpha}\right]~{}~{}\text{for}~{}\alpha\in(0,1)~{}\text{and}~{}\forall~{}t\in[a,b],$

where $E_{\alpha}(\cdot)$ is the one parameter Mittag-Leffler function [30,
Section 1.2].

######  Lemma 2.3.

[7] Consider the following initial value problem

$\displaystyle\partial^{\alpha}_{t}y(t)$
$\displaystyle=g(t,y(t)),~{}t\in(0,T],~{}\alpha\in(0,1),$ (2.19)
$\displaystyle y(0)$ $\displaystyle=y_{0}.$

Let $y_{0}\in\mathbb{R},K^{\ast}>0,t^{\ast}>0$. Define
$D=\left\\{(t,y(t));~{}t\in[0,t^{\ast}],~{}|y-y_{0}|\leq K^{\ast}\right\\}$.
Let function $g:D\rightarrow\mathbb{R}$ be a continuous. Define
$M^{\ast}=\sup_{(t,y(t))\in D}|g(t,y(t)|.$ Then there exists a continuous
function $y\in C[0,T^{\ast}]$ which solves the problem (2.19), where

$T^{\ast}=\begin{cases}t^{\ast};&~{}~{}~{}~{}M^{\ast}=0,\\\
\min\\{t^{\ast},\left(\frac{K^{\ast}\Gamma(1+\alpha)}{M^{\ast}}\right)^{\frac{1}{\alpha}}\\};&~{}~{}~{}~{}\text{else}.\end{cases}$
(2.20)

######  Lemma 2.4.

[18, Lemma 4.1] For $T>0$ and $\alpha\in(0,1)$. Let $X,Y,$ and $Z$ be the
Banach spaces such that $X$ is compactly embedded in $Y$ and $Y$ is
continuously embedded in $Z$. Suppose that $W\subset
L^{1}_{\text{loc}}(0,T;X)$ satisfies the following

  1. 1.

There exist a constant $C_{1}>0$ such that for all $u\in W$

$\sup_{t\in(0,T)}\left(\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}\|u(s)\|^{2}_{X}~{}ds\right)\leq
C_{1}.$ (2.21)

  2. 2.

There exists a constant $C_{2}>0$ such that for all $u\in W$

$\|\partial^{\alpha}_{t}u\|_{L^{2}(0,T;Z)}\leq C_{2}.$ (2.22)

Then $W$ is relatively compact in $L^{2}(0,T;Y)$.

######  Lemma 2.5.

[34] Let $k$ be the kernel defined in (2.16) then there exists a sequence of
kernels $k_{n}$ in $W^{1,1}(0,T)$ such that $k_{n}$ is nonnegative and
nonincreasing in $(0,\infty)$. Also

$k_{n}\rightarrow
k~{}\text{in}~{}L^{1}(0,T)~{}\text{as}~{}n\rightarrow\infty,$ (2.23)

and

$\frac{d}{dt}\left(k_{n}\ast u\right)\rightarrow\frac{d}{dt}(k\ast
u)~{}\text{in}~{}L^{2}(0,T;L^{2}(\Omega))~{}\text{as}~{}n\rightarrow\infty.$
(2.24)

##  3 Main results

The connection between Caputo fractional derivative and Riemann-Liouville
fractional derivative reduces the problem ($D_{\alpha}$) into

$\displaystyle\partial^{\alpha}_{t}u-\nabla\cdot\left(M\left(x,t,\|\nabla
u\|^{2}\right)\nabla u\right)$
$\displaystyle=f(x,t)+\int_{0}^{t}b(x,t,s)u(s)~{}ds~{}\text{in}~{}\Omega\times(0,T],$
($R_{\alpha}$)

with initial and boundary conditions

$\displaystyle u(x,0)$ $\displaystyle=u_{0}(x)\quad\text{in}~{}\Omega,$
$\displaystyle u(x,t)$
$\displaystyle=0\quad\text{on}~{}\partial\Omega\times[0,T].$

The notion of weak formulation for the problem ($R_{\alpha}$) is described as
follows:
find $u$ in $L^{\infty}\left(0,T;L^{2}(\Omega)\right)\cap
L^{2}\left(0,T;H^{1}_{0}(\Omega)\right)$ and $\partial^{\alpha}_{t}u$ in
$L^{2}\left(0,T;{L^{2}(\Omega)}\right)$ such that the following equations hold
for all $v$ in $H^{1}_{0}(\Omega)$ and $a.e.~{}t$ in $(0,T]$

$\displaystyle\left(\partial^{\alpha}_{t}u,v\right)+\left(M\left(x,t,\|\nabla
u\|^{2}\right)\nabla u,\nabla v\right)$
$\displaystyle=\left(f,v\right)+\int_{0}^{t}B(t,s,u(s),v)~{}ds,~{}\text{in}~{}\Omega\times(0,T],$
($W_{\alpha}$) $\displaystyle u(x,0)$
$\displaystyle=u_{0}(x)\quad\text{in}~{}\Omega.$

######  Theorem 3.1.

(Well-posedness of the weak formulation ($W_{\alpha}$)) Under the hypotheses
(H1), (H2), and (H3) the problem ($W_{\alpha}$) admits a unique solution that
satisfies the following a priori bounds

$\|u\|_{L^{\infty}\left(0,T;L^{2}(\Omega)\right)}+\|u\|_{L^{2}_{\alpha}\left(0,T;H^{1}_{0}(\Omega)\right)}\lesssim\left(\|\nabla
u_{0}\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}\right),$ (3.1)
$\|u\|_{L^{\infty}(0,T;H^{1}_{0}(\Omega))}+\|u\|_{L^{2}_{\alpha}\left(0,T;H^{2}(\Omega)\right)}\lesssim\left(\|\nabla
u_{0}\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}\right).$ (3.2)

For the semi discrete formulation of the problem ($D_{\alpha}$), we discretize
the domain in the space variable by a conforming FEM [32] and keep the time
direction continuous. Let $\mathbb{T}_{h}$ be a shape regular (non
overlapping), quasi-uniform triangulation of the domain $\Omega$ and $h$ be
the discretization parameter in the space direction. We define a finite
dimensional subspace $X_{h}$ of $H^{1}_{0}(\Omega)$ as

$X_{h}:=\\{v_{h}\in C(\bar{\Omega})~{}:~{}v_{h}|_{\tau}~{}\text{is a linear
polynomial for
all}~{}\tau\in\mathbb{T}_{h}~{}\text{and}~{}v_{h}=0~{}\text{on}~{}\partial\Omega\\}.$

The semi discrete formulation for the problem ($D_{\alpha}$) is to seek
$u_{h}$ in $X_{h}$ such that the following equations hold for all $v_{h}$ in
$X_{h}$ and $a.e.$ $t$ in $(0,T]$

$\displaystyle\left(\partial^{\alpha}_{t}u_{h},v_{h}\right)$
$\displaystyle+\left(M\left(x,t,\|\nabla u_{h}\|^{2}\right)\nabla u_{h},\nabla
v_{h}\right)$ ($S_{\alpha}$)
$\displaystyle=(f,v_{h})+\int_{0}^{t}B(t,s,u_{h}(s),v_{h})~{}ds~{}\text{in}~{}\mathbb{T}_{h}\times(0,T],$
$\displaystyle u_{h}(x,0)$
$\displaystyle=u_{h}^{0}~{}\text{in}~{}\mathbb{T}_{h},$

where initial condition $u_{h}^{0}$ is in $X_{h}$ which will be chosen later
in the proof of Theorem 3.2.

######  Theorem 3.2.

(Error estimate for the semi discrete formulation ($S_{\alpha}$)) Suppose that
hypotheses (H1), (H2), and (H3) hold. Then we have the following error
estimate for the solution $u_{h}$ of the semi discrete scheme ($S_{\alpha}$)

$\|u-u_{h}\|_{L^{\infty}\left(0,T;L^{2}(\Omega)\right)}+\|u-u_{h}\|_{L^{2}_{\alpha}\left(0,T;H^{1}_{0}(\Omega)\right)}\lesssim
h,$ (3.3)

provided that $u(t)$ is in $H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ for a.e. $t$
in $[0,T]$.

Further, we move to the fully discrete formulation of the problem
($D_{\alpha}$) for that we divide the interval $[0,T]$ into sub intervals of
uniform step size $k$ and $t_{n}=nk$ for $n=0,1,2,3,\dots,N$ with $t_{N}=T$.
We approximate the Caputo fractional derivative by L1 scheme, Kirchhoff type
nonlinearity by linearization, and memory term by modified Simpson’s rule as
follows
L1 approximation scheme [17]: In this scheme, Caputo fractional derivative is
approximated at the point $t_{n}$ using linear interpolation or backward Euler
difference formula as follows

$\displaystyle{}^{C}D^{\alpha}_{t_{n}}u{}$
$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t_{n}}\frac{1}{(t_{n}-s)^{\alpha}}\frac{\partial
u}{\partial s}ds$ (3.4)
$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\sum_{j=1}^{n}\frac{u^{j}-u^{j-1}}{k}\int_{t_{j-1}}^{t_{j}}\frac{1}{(t_{n}-s)^{\alpha}}~{}ds+\mathbb{Q}^{n}$
$\displaystyle=\frac{k^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n}a_{n-j}\left(u^{j}-u^{j-1}\right)+\mathbb{Q}^{n}$
$\displaystyle=\mathbb{D}^{\alpha}_{t}u^{n}+\mathbb{Q}^{n}$

where $a_{i}=(i+1)^{1-\alpha}-i^{1-\alpha},~{}i\geq 0$, $u^{j}=u(x,t_{j})$,
and $\mathbb{Q}^{n}$ is the truncation error.
Linearization: For nonlinear term we use the following linearized
approximation of $u$ at $t_{n}$ given by

$\displaystyle u^{n}$ $\displaystyle\approx
2u^{n-1}-u^{n-2},~{}\text{for}~{}n~{}\geq~{}2$ (3.5)
$\displaystyle:=\bar{u}^{n-1}.$

Modified Simpson’s rule [29]: Let $m_{1}=[k^{-1/2}]$, where $[\cdot]$ denotes
the greatest integer function. Set $k_{1}=m_{1}k$ and $\bar{t}_{j}=jk_{1}$.
Let $j_{n}$ be the largest even integer such that $\bar{t}_{j_{n}}<t_{n}$ and
introduce quadrature points

$\bar{t}_{j}^{n}=\begin{cases}&jk_{1},\quad 0\leq j\leq j_{n},\\\
&\bar{t}_{j}^{n}+(j-j_{n})k,\quad j_{n}\leq j\leq J_{n},\end{cases}$

where $\bar{t}_{J_{n}}^{n}=t_{n-1}$. Then quadrature rule for any function $g$
is as follows

$\displaystyle\int_{0}^{t_{n}}g(s)~{}ds$
$\displaystyle=\sum_{j=0}^{n-1}w_{nj}g(t_{j})+q^{n}(g)$ (3.6)
$\displaystyle=\frac{k_{1}}{3}\sum_{j=1}^{j_{n}/2}\left[g(\bar{t}_{2j}^{n})+4g(\bar{t}_{2j-1}^{n})+g(\bar{t}_{2j-2}^{n})\right]$
$\displaystyle+\frac{k}{2}\sum_{j=j_{n}+1}^{J_{n}}\left[g(\bar{t}_{j}^{n})+g(\bar{t}_{j-1}^{n})\right]+kg(\bar{t}_{J_{n}}^{n})+q^{n}(g),$

where $w_{nj}$ are called quadrature weights and $q^{n}(g)$ is the quadrature
error associated with the function $g$ at $t_{n}$.

On the basis of approximations (3.4), (3.5), and (3.6) we develop the
following linearized L1 Galerkin FEM.
Linearized L1 Galerkin FEM: Find $u_{h}^{n}~{}(n=1,2,3,\dots,N)$ in $X_{h}$
with $\bar{u}_{h}^{n-1}=2u_{h}^{n-1}-u_{h}^{n-2}$ such that the following
equations hold for all $v_{h}$ in $X_{h}$
For $n\geq 2,$

$\left(\mathbb{D}^{\alpha}_{t}u_{h}^{n},v_{h}\right)+\left(M\left(x,t_{n},\|\nabla\bar{u}_{h}^{n-1}\|^{2}\right)\nabla
u_{h}^{n},\nabla
v_{h}\right)=\left(f^{n},v_{h}\right)+\sum_{j=1}^{n-1}w_{nj}B\left(t_{n},t_{j},u_{h}^{j},v_{h}\right).$
($E_{\alpha}$)

For $n=1$,

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}u_{h}^{1},v_{h}\right)+\left(M\left(x,t_{1},\|\nabla
u_{h}^{1}\|^{2}\right)\nabla u_{h}^{1},\nabla
v_{h}\right)=\left(f^{1},v_{h}\right)+kB\left(t_{1},t_{0},u_{h}^{0},v_{h}\right),$

with initial condition $u_{h}^{0}$ that is to be chosen later in the proof of
Theorem 3.2.
To access the convergence rate of the developed numerical scheme
($E_{\alpha}$), we need the following discrete kernel corresponding to the
kernel $(a_{j})$

######  Lemma 3.3.

[17] Let $p_{n}$ be a sequence defined by

$p_{0}=1,~{}~{}p_{n}=\sum_{j=1}^{n}(a_{j-1}-a_{j})p_{n-j}~{}~{}\text{for}~{}~{}n\geq
1.$

Then $p_{n}$ satisfies

$0<p_{n}<1,$ (3.7) $\sum_{j=k}^{n}p_{n-j}a_{j-k}=1,\quad 1\leq k\leq n,$ (3.8)
$\Gamma(2-\alpha)\sum_{j=1}^{n}p_{n-j}\leq\frac{n^{\alpha}}{\Gamma(1+\alpha)}.$
(3.9)

######  Theorem 3.4.

(Convergence estimate for the numerical scheme ($E_{\alpha}$)) Under the
hypotheses (H1), (H2), and (H3) the fully discrete solution
$u_{h}^{n}~{}(1\leq n\leq N)$ of the scheme ($E_{\alpha}$) converges to the
solution $u$ of the problem ($D_{\alpha}$) with the following rate of accuracy

$\max_{1\leq n\leq
N}\|u(t_{n})-u_{h}^{n}\|+\left(k^{\alpha}\sum_{n=1}^{N}p_{N-n}\|\nabla
u(t_{n})-\nabla u_{h}^{n}\|^{2}\right)^{1/2}\lesssim(h+k^{2-\alpha}).$ (3.10)

At this point one can see that convergence rate is of $O(k^{2-\alpha})$ in the
temporal direction. To improve this convergence rate a new linearized
fractional Crank-Nicolson-Galerkin FEM is proposed. In this scheme we replace
the L1 approximation of the Caputo fractional derivative with L2-1σ
$(\sigma=\frac{\alpha}{2})$ scheme [1] at $t_{n-\sigma}$
$(t_{n-\sigma}=(1-\sigma)t_{n}+\sigma t_{n-1})$, linearization technique for
nonlinearity at $t_{n-\sigma}$, and modified Simpson’s rule for the memory
term at $t_{n-\sigma}$.
L2-1σ approximation scheme [1]: In this scheme, Caputo fractional derivative
is approximated at the point $t_{n-\sigma}$ as follows

$^{C}D^{\alpha}_{t_{n-\sigma}}=\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}u+\tilde{\mathbb{Q}}^{n-\sigma},$
(3.11)

where

$\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}u=\frac{k^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n}\tilde{c}_{n-j}^{(n)}\left(u^{j}-u^{j-1}\right),$
(3.12)

with weights $\tilde{c}_{n-j}^{(n)}$ satisfying
$\tilde{c}_{0}^{(1)}=\tilde{a}_{0}$ for $n=1$ and for $n\geq 2$

$\tilde{c}_{j}^{(n)}=\begin{cases}\tilde{a}_{0}+\tilde{b}_{1},&j=0,\\\
\tilde{a}_{j}+\tilde{b}_{j+1}-\tilde{b}_{j},&1\leq j\leq n-2,\\\
\tilde{a}_{j}-\tilde{b}_{j},&j=n-1,\end{cases}$ (3.13)

where

$\tilde{a}_{0}=(1-\sigma)^{1-\alpha}~{}\text{and}~{}\tilde{a}_{l}=\left(l+1-\sigma\right)^{1-\alpha}-\left(1-\sigma\right)^{1-\alpha}\quad
l\geq 1,$
$\tilde{b}_{l}=\frac{1}{(2-\alpha)}\left[\left(l+1-\sigma\right)^{2-\alpha}-\left(l-\sigma\right)^{2-\alpha}\right]-\frac{1}{2}\left[\left(l+1-\sigma\right)^{1-\alpha}+\left(l-\sigma\right)^{1-\alpha}\right]~{}l\geq
1,$

with $\tilde{\mathbb{Q}}^{n-\sigma}$ is the truncation error.
Linearization: Linearized approximation of the nonlinearity and diffusion at
$t_{n-\sigma}$ is given below. For Kirchhoff term

$\displaystyle u^{n-\sigma}$
$\displaystyle\approx(2-\sigma)u^{n-1}-(1-\sigma)u^{n-2},~{}\text{for}~{}n~{}\geq~{}2$
(3.14) $\displaystyle:=\bar{u}^{n-1,\sigma},$

and for diffusion term

$\displaystyle u^{n-\sigma}$
$\displaystyle\approx(1-\sigma)u^{n}+(\sigma)u^{n-1},~{}\text{for}~{}n~{}\geq~{}1$
(3.15) $\displaystyle:=\hat{u}^{n,\sigma}.$

Modified Simpson’s rule: With a small modification in (3.6) we obtain the
following approximation of memory term on $[0,t_{n-\sigma}]$

$\displaystyle\int_{0}^{t_{n-\sigma}}g(s)~{}ds$
$\displaystyle=\sum_{j=0}^{n-1}\tilde{w}_{nj}g(t_{j})+\tilde{q}^{n-\sigma}(g)$
(3.16)
$\displaystyle=\frac{k_{1}}{3}\sum_{j=1}^{j_{n}/2}\left[g(\bar{t}_{2j}^{n})+4g(\bar{t}_{2j-1}^{n})+g(\bar{t}_{2j-2}^{n})\right]$
$\displaystyle+\frac{k}{2}\sum_{j=j_{n}+1}^{J_{n}}\left[g(\bar{t}_{j}^{n})+g(\bar{t}_{j-1}^{n})\right]+\left(1-\sigma\right)kg(\bar{t}_{J_{n}}^{n})+\tilde{q}^{n-\sigma}(g),$

where $\tilde{q}^{n-\sigma}(g)$ is the quadrature error associated with the
function $g$ at $t_{n-\sigma}$.
By combining all approximations (3.11)-(3.16), we construct the following
linearized L2-1σ Galerkin FEM.
Linearized L2-1σ Galerkin FEM: Find $u_{h}^{n}~{}(n=1,2,3,\dots,N)$ in $X_{h}$
with
$\bar{u}_{h}^{n-1,\sigma}=\left(2-\sigma\right)u_{h}^{n-1}-\left(1-\sigma\right)u_{h}^{n-2}$
and
$\hat{u}_{h}^{n,\sigma}=\left(1-\sigma\right)u_{h}^{n}+\left(\sigma\right)u_{h}^{n-1}$
such that the following equations hold for all $v_{h}$ in $X_{h}$
For $n\geq 2$,

$\displaystyle\left(\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}u_{h}^{n},v_{h}\right)+\left(M\left(x,t_{n-\sigma},\|\nabla\bar{u}_{h}^{n-1,\sigma}\|^{2}\right)\nabla\hat{u}_{h}^{n,\sigma},\nabla
v_{h}\right)$
$\displaystyle=\sum_{j=1}^{n-1}\tilde{w}_{nj}B\left(t_{n-\sigma},t_{j},u_{h}^{j},v_{h}\right)$
($F_{\alpha}$) $\displaystyle+\left(f^{n-\sigma},v_{h}\right).$

For $n=1$,

$\displaystyle\left(\tilde{\mathbb{D}}^{\alpha}_{t_{1-\sigma}}u_{h}^{1},v_{h}\right)+\left(M\left(x,t_{1-\sigma},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2}\right)\nabla\hat{u}_{h}^{1,\sigma},\nabla
v_{h}\right)$
$\displaystyle=\left(1-\sigma\right)kB\left(t_{1-\sigma},t_{0},u_{h}^{0},v_{h}\right)$
$\displaystyle+\left(f^{1-\sigma},v_{h}\right),$

Similar to the Lemma 3.3 we have the following discrete kernel corresponding
to the kernel $(\tilde{c}^{n}_{j})$.

######  Lemma 3.5.

[19] Define

$\tilde{p}_{0}^{(n)}=\frac{1}{\tilde{c}_{0}^{(n)}},~{}~{}\tilde{p}_{j}^{(n)}=\frac{1}{\tilde{c}_{0}^{(n-j)}}\sum_{k=0}^{j-1}\left(\tilde{c}_{j-k-1}^{(n-k)}-\tilde{c}_{j-k}^{(n-k)}\right)\tilde{p}_{k}^{(n)}~{}~{}\text{for}~{}~{}1\leq
j\leq n-1.$

Then $\tilde{p}_{j}^{(n)}$ satisfies

$0<\tilde{p}_{n-j}^{(n)}<1,$ (3.17)
$\sum_{j=k}^{n}\tilde{p}_{n-j}^{(n)}\tilde{c}_{j-k}^{(j)}=1,\quad 1\leq k\leq
n\leq N,$ (3.18)
$\Gamma(2-\alpha)\sum_{j=1}^{n}\tilde{p}_{n-j}^{(n)}\leq\frac{n^{\alpha}}{\Gamma(1+\alpha)},\quad
1\leq n\leq N.$ (3.19)

######  Theorem 3.6.

(Convergence estimate for the numerical scheme ($F_{\alpha}$)) Suppose that
hypotheses (H1), (H2), and (H3) hold. Then the fully discrete solution
$u_{h}^{n}~{}(1\leq n\leq N)$ of the scheme ($F_{\alpha}$) satisfies the
following convergence estimate

$\max_{1\leq n\leq
N}\|u(t_{n})-u_{h}^{n}\|+\left(k^{\alpha}\sum_{n=1}^{N}\tilde{p}_{N-n}^{(N)}\|\nabla
u(t_{n})-\nabla u_{h}^{n}\|^{2}\right)^{1/2}\lesssim(h+k^{2}).$ (3.20)

##  4 Well-posedness of the weak formulation ($W_{\alpha}$)

In this section, we prove the well-posedness of the weak formulation
($W_{\alpha}$) using the Galerkin method. For this, we study the variational
formulation ($W_{\alpha}$) and apply the energy argument to derive a priori
bounds on every Galerkin sequence. As a consequence of compactness Lemma 2.4,
these a priori bounds establish the convergence of the Galerkin sequence to
the weak solution of the problem ($D_{\alpha}$).

###  4.1 Proof of the Theorem 3.1

###### Proof.

Let $\\{\lambda_{i}\\}_{i=1}^{m}$ be the eigenvalues and
$\\{\phi_{i}\\}_{i=1}^{m}$ be the corresponding eigenfunctions of the
Dirichlet problem for the standard Laplacian operator in $H^{1}_{0}(\Omega)$.
We consider the finite dimensional subspace $\mathbb{V}_{m}$ of
$H^{1}_{0}(\Omega)$ such that
$\mathbb{V}_{m}=\text{span}\\{\phi_{1},\phi_{2},\dots,\phi_{m}\\}.$ We assume
that

$u_{m}(\cdot,t)=\sum_{j=1}^{m}\alpha_{mj}(t)\phi_{j}\quad\text{and }\quad
u_{m}(\cdot,0)=\sum_{j=1}^{m}(u_{0},\phi_{j})\phi_{j},$ (4.1)

then $u_{m}(\cdot,0)$ converges to $u_{0}$. Consider the problem
($W_{\alpha}$) onto finite dimensional subspace $\mathbb{V}_{m}$, $i.e.,$ find
$u_{m}\in\mathbb{V}_{m}$ which satisfying the following equation for all
$v_{m}\in\mathbb{V}_{m}~{}\text{and}~{}a.e.~{}t\in(0,T]$

$\displaystyle\left(\partial^{\alpha}_{t}u_{m},v_{m}\right)+\left(M(x,t,\|\nabla
u_{m}\|^{2})\nabla u_{m},\nabla
v_{m}\right)=(f,v_{m})+\int_{0}^{t}B(t,s,u_{m}(s),v_{m})~{}ds.$ (4.2)

Put the values of $u_{m}$ and $u_{m}(0)$ in (4.2), we obtain a system of
fractional order differential equations. Then by the theory of fractional
order differential equations Lemma 2.3, the system (4.2) has a continuous
solution $u_{m}(t)$ on some interval $[0,t_{n}),0<t_{n}<T,$ with vanishing
trace of $k\ast(u_{m}-u_{m}(0))$ at $t=0$ [34]. These local solutions
$u_{m}(t)$ are extended to the whole interval by using the following a priori
bounds.
(A priori bounds) Take $v_{m}=u_{m}(t)$ in (4.2) to get

$\displaystyle\left(\frac{d}{dt}\left(k\ast u_{m}\right),u_{m}\right)$
$\displaystyle+\left(M\left(x,t,\|\nabla u_{m}\|^{2}\right)\nabla u_{m},\nabla
u_{m}\right)$ (4.3) $\displaystyle=\left(\frac{d}{dt}\left(k\ast
u_{m}(0)\right),u_{m}\right)+(f,u_{m})+\int_{0}^{t}B(t,s,u_{m}(s),u_{m}(t))~{}ds.$

Let $k_{n},n\in\mathbb{N}$ be the sequence of kernels defined in Lemma 2.5,
then equation (4.3) is rewritten as

$\displaystyle\left(\frac{d}{dt}(k_{n}\ast
u_{m}),u_{m}\right)+\left(M\left(x,t,\|\nabla u_{m}\|^{2}\right)\nabla
u_{m},\nabla u_{m}\right)$ (4.4)
$\displaystyle=h_{mn}(t)+k(t)(u_{m}(0),u_{m})+(f,u_{m})+\int_{0}^{t}B(t,s,u_{m}(s),u_{m}(t))~{}ds,$

with

$h_{mn}(t)=\left(\frac{d}{dt}\left(k_{n}\ast
u_{m}\right)-\frac{d}{dt}\left(k\ast u_{m}\right),u_{m}\right).$ (4.5)

Use Lemma 2.1, positivity of diffusion coefficient (H2), continuity of
$B(t,s,\cdot,\cdot)$ (2.15), Cauchy-Schwarz and Young’s inequalities to obtain

$\displaystyle\frac{d}{dt}\left(k_{n}\ast\|u_{m}\|^{2}\right)(t)+k_{n}(t)\|u_{m}\|^{2}+\|\nabla
u_{m}\|^{2}$ (4.6) $\displaystyle\lesssim
h_{mn}(t)+k(t)\|u_{m}(0)\|^{2}+k(t)\|u_{m}\|^{2}+\|f\|^{2}+\|u_{m}\|^{2}+\int_{0}^{t}\|\nabla
u_{m}(s)\|^{2}~{}ds.$

By convolving the equation (4.6) with the kernel
$l(t)=\frac{t^{\alpha-1}}{\Gamma(\alpha)}$ and letting $n\rightarrow\infty$ in
(4.6) then equation (4.6) reduces to

$\displaystyle\|u_{m}(t)\|^{2}+\left(l\ast\|\nabla u_{m}\|^{2}\right)(t)$
$\displaystyle\lesssim\|u_{m}(0)\|^{2}+\left(l\ast\|f\|^{2}\right)(t)$ (4.7)
$\displaystyle+l\ast\left(\|u_{m}\|^{2}+\int_{0}^{t}\|\nabla
u_{m}(s)\|^{2}~{}ds\right).$

In (4.7) we have used the fact that

$l\ast\frac{d}{dt}\left(k_{n}\ast\|u_{m}\|^{2}\right)(t)=\frac{d}{dt}\left(k_{n}\ast
l\ast\|u_{m}\|^{2}\right)(t)\rightarrow\frac{d}{dt}\left(k\ast
l\ast\|u_{m}\|^{2}\right)(t)=\|u_{m}\|^{2},$ (4.8)

with

$(l\ast k_{n})(t)\rightarrow(l\ast k)(t)=1~{}~{}\text{and}~{}~{}(l\ast
h_{mn})(t)\rightarrow
0~{}~{}\text{as}~{}~{}n\rightarrow\infty~{}~{}\text{in}~{}~{}L^{1}(0,T).$
(4.9)

Denote $\tilde{u}_{m}(t):=\|u_{m}(t)\|^{2}+\left(l\ast\|\nabla
u_{m}\|^{2}\right)(t)$ and
$\tilde{v}_{m}(t):=\|u_{m}(0)\|^{2}+\left(l\ast\|f\|^{2}\right)(t)$. Then
equation (4.7) is converted into

$\displaystyle\tilde{u}_{m}(t)$
$\displaystyle\lesssim\tilde{v}_{m}(t)+l\ast\left(\|u_{m}\|^{2}+\int_{0}^{t}\|\nabla
u_{m}(s)\|^{2}~{}ds\right)$ (4.10)
$\displaystyle\lesssim\tilde{v}_{m}(t)+l\ast\left(\|u_{m}\|^{2}+\int_{0}^{t}(t-s)^{\alpha-1}(t-s)^{1-\alpha}\|\nabla
u_{m}(s)\|^{2}~{}ds\right)$
$\displaystyle\lesssim\tilde{v}_{m}(t)+\int_{0}^{t}(t-s)^{\alpha-1}\tilde{u}_{m}(s)~{}ds.$

As a consequence of Lemma 2.2 and Poincaré inequality, we deduce

$\displaystyle\|u_{m}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|u_{m}\|_{L^{2}_{\alpha}(0,T;H^{1}_{0}(\Omega))}$
$\displaystyle\lesssim\|u_{m}(0)\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}$
(4.11) $\displaystyle\lesssim\|u_{0}\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}$
$\displaystyle\lesssim\|\nabla u_{0}\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}.$

Further, we define two new discrete Laplacian operators
$\Delta_{m}^{M},~{}\Delta_{m}^{b_{2}}:\mathbb{V}_{m}\rightarrow\mathbb{V}_{m}$
such that

$(-\Delta_{m}^{M}~{}u_{m},v_{m}):=(M(x,t,\|\nabla u_{m}\|^{2})\nabla
u_{m},\nabla v_{m})~{}\text{for
all}~{}u_{m},v_{m}~{}\in~{}\mathbb{V}_{m},~{}t\in(0,T],$ (4.12)

and

$(-\Delta_{m}^{b_{2}}~{}u_{m},v_{m}):=(b_{2}(x,t,s)\nabla u_{m},\nabla
v_{m})~{}\text{for
all}~{}u_{m},v_{m}~{}\in~{}\mathbb{V}_{m},~{}t,s~{}\in(0,T].$ (4.13)

Since diffusion coefficient is positive and $b_{2}$ is a symmetric positive
definite matrix therefore $\Delta_{m}^{M}$ and $\Delta_{m}^{b_{2}}$ are well
defined. We make use of these definitions (4.12) and (4.13) to convert the
equation (4.2) into

$\displaystyle\left(\partial^{\alpha}_{t}u_{m},v_{m}\right)+\left(-\Delta_{m}^{M}~{}u_{m},v_{m}\right)$
(4.14)
$\displaystyle=\left(f,v_{m}\right)+\int_{0}^{t}(-\Delta_{m}^{b_{2}}~{}u_{m}(s),v_{m})~{}ds+\int_{0}^{t}(\nabla\cdot(b_{1}(x,t,s)u_{m}(s)),v_{m})~{}ds$
$\displaystyle+\int_{0}^{t}(b_{0}(x,t,s)u_{m}(s),v_{m})~{}ds.$

Put $v_{m}=-\Delta_{m}^{M}~{}u_{m}(t)$ in (4.14) and apply Cauchy-Schwarz
inequality together with Young’s inequality to obtain

$\displaystyle\left(\partial^{\alpha}_{t}\nabla u_{m},M\left(x,t,\|\nabla
u_{m}\|^{2}\right)\nabla u_{m}\right)$
$\displaystyle+\|\Delta_{m}^{M}~{}u_{m}\|^{2}$ (4.15)
$\displaystyle\lesssim\int_{0}^{t}\|\Delta_{m}^{b_{2}}~{}u_{m}(s)\|^{2}~{}ds+\int_{0}^{t}\|\nabla
u_{m}(s)\|^{2}~{}ds$
$\displaystyle+\int_{0}^{t}\|u_{m}(s)\|^{2}~{}ds+\|f\|^{2}.$

Estimate $|(b_{2}(x,t,s)\nabla u_{m}(s),\nabla v_{m})|\lesssim\|\nabla
u_{m}(s)\|~{}\|\nabla v_{m}\|$ implies

$\displaystyle\|\Delta_{m}^{b_{2}}~{}u_{m}(s)\|=\sup_{v_{m}\in\mathbb{V}_{m}}\frac{|(b_{2}(x,t,s)\nabla
u_{m}(s),\nabla v_{m})|}{\|\nabla v_{m}\|}$ $\displaystyle\lesssim\|\nabla
u_{m}(s)\|.$ (4.16)

Hypothesis (H2) and estimate (4.16) yield

$\displaystyle\left(\frac{d}{dt}\left[k\ast\nabla u_{m}\right](t),\nabla
u_{m}\right)+\|\Delta_{m}^{M}~{}u_{m}\|^{2}$
$\displaystyle\lesssim\int_{0}^{t}\left(\|\nabla
u_{m}(s)\|^{2}+\|u_{m}(s)\|^{2}\right)~{}ds+\|f\|^{2}$ (4.17)
$\displaystyle+k(t)\left(\nabla u_{m}(0),\nabla u_{m}\right).$

Following the similar lines of the proof of estimate (4.11), we reach at

$\displaystyle\|u\|^{2}_{L^{\infty}(0,T;H^{1}_{0}(\Omega))}+\left(l\ast\|\Delta_{m}^{M}~{}u_{m}\|^{2}\right)(t)$
$\displaystyle\lesssim\|\nabla
u_{m}(0)\|^{2}+\|f\|^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}$ (4.18)
$\displaystyle\lesssim\|\nabla
u_{0}\|^{2}+\|f\|^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}.$

Finally, substitute $v_{m}=\partial^{\alpha}_{t}u_{m}$ in (4.14) to have

$\displaystyle\|\partial^{\alpha}_{t}u_{m}\|^{2}+\left(-\Delta_{m}^{M}~{}u_{m},\partial^{\alpha}_{t}u_{m}\right)$
$\displaystyle\lesssim\|f\|^{2}+\int_{0}^{t}(\|\Delta_{m}^{b_{2}}~{}u_{m}(s)\|^{2}+\|\nabla
u_{m}(s)\|^{2})~{}ds.$ (4.19)
$\displaystyle+\int_{0}^{t}\|u_{m}(s)\|^{2}~{}ds.$

Proceeding further as estimate (4.18) is proved to conclude

$\displaystyle\|\partial^{\alpha}_{t}u_{m}\|_{L^{2}(0,T;L^{2}(\Omega))}+\|u_{m}\|_{L^{\infty}(0,T;H^{1}_{0}(\Omega))}$
$\displaystyle\lesssim\|\nabla
u_{m}(0)\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ (4.20)
$\displaystyle\lesssim\|\nabla u_{0}\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}.$

Thus, estimates (4.11) and (4.20) provide a subsequence of $(u_{m})$ again
denoted by $(u_{m})$ such that $u_{m}\rightharpoonup u$ in
$L^{2}\left(0,T;H^{1}_{0}(\Omega)\right)$ and
$\partial^{\alpha}_{t}u_{m}\rightharpoonup\partial^{\alpha}_{t}u$ in
$L^{2}\left(0,T;L^{2}(\Omega)\right)$. In the light of estimates (4.18) and
(4.20), we apply compactness Lemma 2.4 to conclude $u_{m}\rightarrow u$ in
$L^{2}\left(0,T;H^{1}_{0}(\Omega)\right)$. Now using the fact that
$M(x,t,\|\nabla u_{m}\|^{2})$ and $B(t,s,u_{m}(s),v_{m})$ are continuous and
an application of Lebesgue dominated convergence theorem, we pass the limit
inside (4.2) which establishes the existence of weak solutions of the problem
($D_{\alpha}$).
(Initial Condition) The weak solution $u$ satisfies the following equation for
all $v$ in $H^{1}_{0}(\Omega)$

$\displaystyle\left(\frac{d}{dt}\left[k\ast(u-u_{0})\right](t),v\right)+\left(M\left(x,t,\|\nabla
u\|^{2}\right)\nabla u,\nabla v\right)$
$\displaystyle=\left(f,v\right)+\int_{0}^{t}B(t,s,u(s),v)ds.$ (4.21)

Let $\phi$ in $C^{1}\left([0,T];H^{1}_{0}(\Omega)\right)$ with $\phi(T)=0$,
multiply (4.21) with $\phi$ and integrate by parts to get

$\displaystyle-\int_{0}^{T}\left((k\ast(u-u_{0}))(t),v\right)\phi^{\prime}(t)dt$
$\displaystyle+\int_{0}^{T}\left(M\left(x,t,\|\nabla u\|^{2}\right)\nabla
u,\nabla v\right)\phi(t)dt$ (4.22)
$\displaystyle=\int_{0}^{T}\left(f,v\right)\phi(t)dt+\int_{0}^{T}\int_{0}^{t}B(t,s,u(s),v)\phi(t)dsdt$
$\displaystyle+((k\ast(u-u_{0}))(0),\phi(0)).$

Since $C^{1}\left([0,T];H^{1}_{0}(\Omega)\right)$ is dense in
$L^{2}\left(0,T;H^{1}_{0}(\Omega)\right)$, thus using (4.2) and
$(k\ast(u_{m}-u_{m}(0)))$ has vanishing trace at $t=0$, we have

$\displaystyle-\int_{0}^{T}\left((k\ast(u_{m}-u_{m}(0)))(t),v\right)$
$\displaystyle\phi^{\prime}(t)dt+\int_{0}^{T}\left(M\left(x,t,\|\nabla
u_{m}\|^{2}\right)\nabla u_{m},\nabla v\right)\phi(t)dt$ (4.23)
$\displaystyle=\int_{0}^{T}\left(f,v\right)\phi(t)dt+\int_{0}^{T}\int_{0}^{t}B(t,s,u_{m}(s),v)\phi(t)dsdt.$

Let $m$ tend to infinity in (4.23) and comparing with (4.22) to get
$((k\ast(u-u_{0}))(0),\phi(0))=0$. Since $\phi(0)$ is arbitrary, so we have
$(k\ast(u-u_{0}))(0)=0$ which implies $u=u_{0}$ at $t=0$ [34].
(Uniqueness) Suppose that $u_{1},u_{2}$ are solutions of the weak formulation
($W_{\alpha}$), then $z=u_{1}-u_{2}$ satisfies the following equation for all
$v\in H_{0}^{1}(\Omega)$ and $a.e.~{}t\in(0,T]$

$\displaystyle\left(\frac{d}{dt}(k\ast z)(t),v\right)$
$\displaystyle+\left(M\big{(}x,t,\|\nabla u_{1}\|^{2}\big{)}\nabla z,\nabla
v\right)$ (4.24) $\displaystyle=\left(\left[M\big{(}x,t,\|\nabla
u_{2}\|^{2}\big{)}-M\big{(}x,t,\|\nabla u_{1}\|^{2}\big{)}\right]\nabla
u_{2},\nabla v\right)$ $\displaystyle+\int_{0}^{t}B(t,s,z(s),v)~{}ds.$

Put $v=z(t)$ in (4.24) and using (H2), (H3), and a priori bound (3.2) on
$u_{1},u_{2}$ along with Cauchy-Schwarz and Young’s inequality to obtain

$\displaystyle\left(\frac{d}{dt}(k\ast
z)(t),z(t)\right)+(m_{0}-4L_{M}K^{2})\|\nabla
z\|^{2}\lesssim\int_{0}^{t}\|\nabla z(s)\|^{2}ds.$ (4.25)

Following the similar lines as in the proof of estimate (4.11) and using (H2)
we conclude
$\|z\|_{L^{2}\left(0,T;H^{1}_{0}(\Omega)\right)}=\|z\|_{L^{\infty}\left(0,T;L^{2}(\Omega)\right)}=0$.
Thus uniqueness follows. ∎

##  5 Semi discrete formulation and error estimate

In this section, we discuss the well-posedness of the semi discrete
formulation ($S_{\alpha}$) and derive error estimate for the semi discrete
solution by modifying Ritz-Volterra projection operator.

######  Theorem 5.1.

Suppose that hypotheses (H1), (H2), and (H3) hold. Then there exists a unique
solution to the problem ($S_{\alpha}$) which satisfies the following a priori
bounds

$\|u_{h}\|_{L^{\infty}\left(0,T;L^{2}(\Omega)\right)}+\|u_{h}\|_{L^{2}_{\alpha}\left(0,T;H^{1}_{0}(\Omega)\right)}\lesssim\left(\|\nabla
u_{0}\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}\right),$ (5.1)
$\|\partial^{\alpha}_{t}u_{h}\|_{L^{2}(0,T;L^{2}(\Omega))}+\|u_{h}\|_{L^{\infty}(0,T;H^{1}_{0}(\Omega))}\lesssim\left(\|\nabla
u_{0}\|+\|f\|_{L^{\infty}(0,T;L^{2}(\Omega))}\right).$ (5.2)

###### Proof.

This theorem is proved analogously to the proof of Theorem 3.1. ∎

For the semi discrete error estimate, we define a new Ritz-Volterra type
projection operator $W:[0,T]\rightarrow X_{h}$ by

$\left(M\left(x,t,\|\nabla u\|^{2}\right)\nabla(u-W),\nabla
v_{h}\right):=\int_{0}^{t}B(t,s,u(s)-W(s),v_{h})~{}ds\quad\text{for
all}~{}v_{h}~{}\text{in}~{}X_{h}.$ (5.3)

This modified Ritz-Volterra projection operator $W$ is well defined by the
positivity of the Kirchhoff term $M$ [5]. This projection operator satisfies
the following stability and best approximation results.

######  Lemma 5.2.

[16] Consider $W$ is the modified Ritz-Volterra projection operator defined in
(5.3), then $\|\nabla W\|$ is bounded for every $t$ in $[0,T]$, i.e.,

$\|\nabla W\|\lesssim\|\nabla u\|.$

To derive the best approximation properties of the modified Ritz-Volterra
projection operator, we assume some regularity assumptions on the solution $u$
of the problem ($D_{\alpha}$) such that

$\|u(t)\|_{2}\lesssim C~{}\text{and}~{}\|u_{t}(t)\|_{2}\lesssim
C~{}\forall~{}t~{}\in~{}[0,T].$ (5.4)

######  Theorem 5.3.

Suppose that the solution $u$ of the problem ($D_{\alpha}$) satisfies (5.4).
Then modified Ritz-Volterra projection operator has the following best
approximation properties

$\displaystyle\|\rho(t)\|+h\|\nabla\rho(t)\|$ $\displaystyle\lesssim
h^{2}~{}\forall~{}t~{}\in~{}[0,T],$ (5.5)
$\displaystyle\|\rho_{t}(t)\|+h\|\nabla\rho_{t}(t)\|$ $\displaystyle\lesssim
h^{2}~{}\forall~{}t~{}\in~{}[0,T],$

where $\rho:=u-W$.

###### Proof.

For the proof of this theorem we refer the readers to [4, 5]. ∎

Now error estimate for the semi discrete formulation ($S_{\alpha}$) is
attained as stated in Theorem 3.2.

###  5.1 Proof of the Theorem 3.2

###### Proof.

Denote $(W-u_{h}):=\theta$ such that $u-u_{h}=\rho+\theta.$ Then put
$u_{h}=W-\theta$ in the problem ($S_{\alpha}$) to have

$\displaystyle\left(\partial^{\alpha}_{t}(W-\theta),v_{h}\right)$
$\displaystyle+\left(M\left(x,t,\|\nabla
u_{h}\|^{2}\right)\nabla(W-\theta),\nabla v_{h}\right)$
$\displaystyle=(f,v_{h})+\int_{0}^{t}B(t,s,W(s)-\theta(s),v_{h})~{}ds.$

Weak formulation ($W_{\alpha}$) and the definition (5.3) of the modified Ritz-
Volterra projection operator $W$ yield

$\displaystyle\left(\partial^{\alpha}_{t}\theta,v_{h}\right)+\left(M(x,t,\|\nabla
u_{h}\|^{2}\right)\nabla\theta,\nabla v_{h})$ (5.6)
$\displaystyle=-\left(\partial^{\alpha}_{t}\rho,v_{h}\right)+\int_{0}^{t}B(t,s,\theta(s),v_{h})ds+\left(\left(M(x,t,\|\nabla
u_{h}\|^{2}-M(x,t,\|\nabla u\|^{2}\right)\nabla W,\nabla v_{h}\right).$

Set $v_{h}=\theta(t)$ in (5.6) and employ (H2), (H3) to obtain

$\displaystyle\left(\partial^{\alpha}_{t}\theta,\theta\right)+m_{0}\|\nabla\theta\|^{2}$
$\displaystyle=\|\partial^{\alpha}_{t}\rho\|\|\theta(t)\|+\|\nabla\theta(t)\|\int_{0}^{t}\|\nabla\theta(s)\|ds$
(5.7) $\displaystyle+L_{M}(\|\nabla u_{h}\|+\|\nabla
u\|)(\|\nabla\rho\|+\|\nabla\theta\|)\|\nabla W\|\|\nabla\theta\|.$

By utilizing proved a priori bounds on $\|\nabla u\|,\|\nabla u_{h}\|$ and
$\|\nabla W\|$ together with Cauchy-Schwarz and Young’s inequality, we obtain

$\displaystyle\left(\partial^{\alpha}_{t}\theta,\theta\right)+(m_{0}-4L_{M}K^{2})\|\nabla\theta\|^{2}$
$\displaystyle\lesssim\|\partial^{\alpha}_{t}\rho\|^{2}+\|\theta\|^{2}+\int_{0}^{t}\|\nabla\theta(s)\|^{2}ds+\|\nabla\rho\|^{2}.$
(5.8)

Use (H2) and apply similar arguments as we prove estimate (4.11) to deduce

$\displaystyle\|\theta\|^{2}_{L^{\infty}\left(0,T;L^{2}(\Omega)\right)}$
$\displaystyle+\left(l\ast\|\nabla\theta\|^{2}\right)(t)\lesssim\left[l\ast\left(\|\nabla\rho\|^{2}+\|\partial^{\alpha}_{t}\rho\|^{2}\right)\right](t)+\|\nabla\theta(0)\|^{2}.$

Consider

$\displaystyle\|\partial^{\alpha}_{t}\rho\|=\|~{}^{C}D^{\alpha}_{t}\rho\|$
$\displaystyle=\left\|\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}\frac{\partial\rho}{\partial
s}(s)~{}ds\right\|$ (5.9)
$\displaystyle\lesssim\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}\left\|\frac{\partial\rho}{\partial
s}(s)\right\|~{}ds$
$\displaystyle\lesssim\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}h^{2}~{}ds\lesssim
h^{2}.$

We choose $u_{h}^{0}=W(0)$ such that $\theta(0)=0$. Then apply approximation
properties of modified Ritz-Volterra projection operator (5.5) and (5.9) to
conclude

$\displaystyle\|\theta\|^{2}_{L^{\infty}\left(0,T;L^{2}(\Omega)\right)}+\left(l\ast\|\nabla\theta\|^{2}\right)(t)\lesssim
h^{2}+h^{4}\lesssim h^{2}.$

Finally, triangle inequality and estimate (5.5) finish the proof. ∎

##  6 Linearized L1 Galerkin FEM

In this section, we prove the well-posedness of the numerical scheme
($E_{\alpha}$) and carry out its convergence analysis. The following two
lemmas provide a priori bounds on the solution of the problem ($E_{\alpha}$).

######  Lemma 6.1.

Under the hypothesis (H1), (H2), and (H3) the solution $u_{h}^{n}~{}(n\geq 1)$
of the scheme ($E_{\alpha}$) satisfy the following a priori bound

$\displaystyle\max_{1\leq m\leq
N}\|u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{N}p_{N-n}\|\nabla
u_{h}^{n}\|^{2}\lesssim\|\nabla u_{0}\|^{2}+\max_{1\leq n\leq
N}\|f^{n}\|^{2}.$ (6.1)

###### Proof.

Put $v_{h}=u_{h}^{1}$ for $n=1$ in the formulation ($E_{\alpha}$) to get

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}u_{h}^{1},u_{h}^{1}\right)+\left(M\left(x,t_{1},\|\nabla
u_{h}^{1}\|^{2}\right)\nabla u_{h}^{1},\nabla
u_{h}^{1}\right)=\left(f^{1},u_{h}^{1}\right)+kB\left(t_{1},t_{0},u_{h}^{0},u_{h}^{1}\right).$

Employing (H2) and (H3) to obtain

$\displaystyle(1-k^{\alpha}\Gamma(2-\alpha))\|u_{h}^{1}\|^{2}+k^{\alpha}\|\nabla
u_{h}^{1}\|^{2}$ $\displaystyle\lesssim
k^{\alpha}\left(\|f^{1}\|^{2}+k^{2}\|\nabla
u_{h}^{0}\|^{2}\right)+\|u_{h}^{0}\|^{2}.$

For sufficiently small $k$ such that $k^{\alpha}<\frac{1}{\Gamma(2-\alpha)}$,
we conclude

$\displaystyle\|u_{h}^{1}\|^{2}+k^{\alpha}\|\nabla u_{h}^{1}\|^{2}$
$\displaystyle\lesssim\|f^{1}\|^{2}+\|\nabla
u_{h}^{0}\|^{2}\lesssim\|f^{1}\|^{2}+\|\nabla u_{0}\|^{2}.$

Further, set $v_{h}=u_{h}^{n}$ for $n\geq 2$ in the scheme ($E_{\alpha}$) to
have

$\left(\mathbb{D}^{\alpha}_{t}u_{h}^{n},u_{h}^{n}\right)+\left(M\left(x,t_{n},\|\nabla\bar{u}_{h}^{n-1}\|^{2}\right)\nabla
u_{h}^{n},\nabla
u_{h}^{n}\right)=\left(f^{n},u_{h}^{n}\right)+\sum_{j=1}^{n-1}w_{nj}B\left(t_{n},t_{j},u_{h}^{j},u_{h}^{n}\right).$
(6.2)

Applying the identity
$\left(\mathbb{D}^{\alpha}_{t}u_{h}^{n},u_{h}^{n}\right)\geq\frac{1}{2}\mathbb{D}^{\alpha}_{t}\|u_{h}^{n}\|^{2}$
[17] and hypotheses (H2), (H3) along with Cauchy-Schwarz and Young’s
inequality to reach at

$\mathbb{D}^{\alpha}_{t}\|u_{h}^{n}\|^{2}+\|\nabla
u_{h}^{n}\|^{2}\lesssim\left(\|f^{n}\|^{2}+\|u_{h}^{n}\|^{2}+\sum_{j=1}^{n-1}w_{nj}\|\nabla
u_{h}^{j}\|^{2}\right).$ (6.3)

By the definition of $\mathbb{D}^{\alpha}_{t}\|u_{h}^{n}\|^{2}$ (3.4), the
equation (6.3) reduces to

$\frac{k^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n}a_{n-j}\left(\|u_{h}^{j}\|^{2}-\|u_{h}^{j-1}\|^{2}\right)+\|\nabla
u_{h}^{n}\|^{2}\lesssim\left(\|f^{n}\|^{2}+\|u_{h}^{n}\|^{2}+\sum_{j=1}^{n-1}w_{nj}\|\nabla
u_{h}^{j}\|^{2}\right).$ (6.4)

Multiply the equation (6.4) by discrete convolution $P_{m-n}$ defined in Lemma
3.3 and take summation from $n=1$ to $m$ to obtain

$\displaystyle\sum_{n=1}^{m}p_{m-n}\sum_{j=1}^{n}a_{n-j}\left(\|u_{h}^{j}\|^{2}-\|u_{h}^{j-1}\|^{2}\right)+k^{\alpha}\Gamma(2-\alpha)\sum_{n=1}^{m}p_{m-n}\|\nabla
u_{h}^{n}\|^{2}$ $\displaystyle\lesssim
k^{\alpha}\Gamma(2-\alpha)\left(\sum_{n=1}^{m}p_{m-n}\|f^{n}\|^{2}+\sum_{n=1}^{m}p_{m-n}\|u_{h}^{n}\|^{2}+\sum_{n=1}^{m}p_{m-n}\sum_{j=1}^{n-1}w_{nj}\|\nabla
u_{h}^{j}\|^{2}\right).$

Interchanging the summation and property of discrete kernel (3.8) with
$\alpha_{0}=k^{\alpha}\Gamma(2-\alpha)$ yield

$\displaystyle(1-\alpha_{0})\|u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|\nabla
u_{h}^{n}\|^{2}$
$\displaystyle\lesssim\sum_{n=1}^{m-1}\left(k^{\alpha}p_{m-n}\|u_{h}^{n}\|^{2}+k_{1}k^{\alpha}\sum_{j=1}^{n}p_{n-j}\|\nabla
u_{h}^{j}\|^{2}\right)$
$\displaystyle+k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|f_{h}^{n}\|^{2}+\|u_{h}^{0}\|^{2}.$

Then for sufficiently small $k^{\alpha}<\frac{1}{\Gamma{(2-\alpha)}}$ one have

$\displaystyle\|u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|\nabla
u_{h}^{n}\|^{2}$
$\displaystyle\lesssim\sum_{n=1}^{m-1}\left(k^{\alpha}p_{m-n}+k_{1}\right)\left(\|u_{h}^{n}\|^{2}+k^{\alpha}\sum_{j=1}^{n}p_{n-j}\|\nabla
u_{h}^{j}\|^{2}\right)$
$\displaystyle+k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|f^{n}\|^{2}+\|u_{h}^{0}\|^{2}.$

Further, the discrete Grönwall’s inequality provides

$\displaystyle\|u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|\nabla
u_{h}^{n}\|^{2}$
$\displaystyle\lesssim\left(k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|f^{n}\|^{2}+\|u_{h}^{0}\|^{2}\right)\exp\left(\sum_{n=1}^{m-1}\left(k^{\alpha}p_{m-n}+k_{1}\right)\right).$
(6.5)

Finally, using property of discrete kernel (3.9) one obtain

$k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|f^{n}\|^{2}\lesssim\max_{1\leq n\leq
N}\|f^{n}\|^{2}\left(k^{\alpha}\sum_{n=1}^{m}p_{m-n}\right)\lesssim\max_{1\leq
n\leq N}\|f^{n}\|^{2}k^{\alpha}m^{\alpha}\lesssim T^{\alpha}\max_{1\leq n\leq
N}\|f^{n}\|^{2}.$ (6.6)

Also,

$\sum_{n=1}^{m-1}\left(k^{\alpha}p_{m-n}+k_{1}\right)\lesssim
k^{\alpha}m^{\alpha}+mk_{1}\lesssim T.$ (6.7)

To conclude the result (6.1), put (6.6) and (6.7) in (6.5) as

$\displaystyle\|u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{m}p_{m-n}\|\nabla
u_{h}^{n}\|^{2}$ $\displaystyle\lesssim\max_{1\leq n\leq
N}\|f^{n}\|^{2}+\|u_{h}^{0}\|^{2}\lesssim\max_{1\leq n\leq
N}\|f^{n}\|^{2}+\|\nabla u_{0}\|^{2}.$ (6.8)

∎

######  Lemma 6.2.

Suppose that hypotheses (H1), (H2), and (H3) hold. Then the solution
$u_{h}^{n}~{}(n\geq 1)$ of the scheme ($E_{\alpha}$) satisfy the following a
priori bound

$\displaystyle\max_{1\leq m\leq N}\|\nabla
u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{N}p_{N-n}\|\Delta_{h}^{M}u_{h}^{n}\|^{2}\lesssim\|\nabla
u_{0}\|^{2}+\max_{1\leq n\leq N}\|f^{n}\|^{2}.$ (6.9)

where $\Delta_{h}^{M}:X_{h}\rightarrow X_{h}$ is the discrete Laplacian
operator defined in (4.12) associated with $M$.

###### Proof.

By making use of definitions of discrete Laplacian operator
$-\Delta_{h}^{M}u_{h}^{1}$ (4.12) and $-\Delta_{h}^{b_{2}}u_{h}^{0}$ (4.13)
the equation for $n=1$ in the scheme ($E_{\alpha}$) is rewritten as

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}u_{h}^{1},v_{h}\right)+(-\Delta_{h}^{M}u_{h}^{1},v_{h})$
$\displaystyle=\left(f^{1},v_{h}\right)+k(-\Delta_{h}^{b_{2}}u_{h}^{0},v_{h})+k(\nabla\cdot(b_{1}(x,t,s)u_{h}^{0}),v_{h})$
(6.10) $\displaystyle+k(b_{0}(x,t,s)u_{h}^{0},v_{h}).$

Setting $v_{h}=-\Delta_{M_{h}}u_{h}^{1}$ in (6.10) one obtain

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}\nabla u_{h}^{1},\nabla
u_{h}^{1}\right)+\|\Delta_{h}^{M}u_{h}^{1}\|^{2}$
$\displaystyle\lesssim\|f^{1}\|^{2}+k^{2}\left(\|\Delta_{h}^{b_{2}}u_{h}^{0}\|^{2}+\|\nabla
u_{h}^{0}\|^{2}+\|u_{h}^{0}\|^{2}\right).$ (6.11)

Identity
$\left(\mathbb{D}^{\alpha}_{t}u_{h}^{n},u_{h}^{n}\right)\geq\frac{1}{2}\mathbb{D}^{\alpha}_{t}\|u_{h}^{n}\|^{2}$
and estimate (4.16) simplify the equation (6.11) to

$\displaystyle\|\nabla
u_{h}^{1}\|^{2}+k^{\alpha}\|\Delta_{h}^{M}u_{h}^{1}\|^{2}$
$\displaystyle\lesssim k^{\alpha}\|f^{1}\|^{2}+k^{\alpha}k^{2}\|\nabla
u_{h}^{0}\|^{2}+\|\nabla u_{h}^{0}\|^{2}.$ (6.12)

For sufficiently small $k$, we deduce

$\displaystyle\|\nabla
u_{h}^{1}\|^{2}+k^{\alpha}\|\Delta_{h}^{M}u_{h}^{1}\|^{2}$
$\displaystyle\lesssim\|f^{1}\|^{2}+\|\nabla
u_{h}^{0}\|^{2}\lesssim\|f^{1}\|^{2}+\|\nabla u_{0}\|^{2}.$ (6.13)

Consider the scheme ($E_{\alpha}$) for $n\geq 2$ with definitions of discrete
Laplacian operators (4.12) and (4.13)

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}u_{h}^{n},v_{h}\right)+\left(-\Delta_{h}^{M}u_{h}^{n},v_{h}\right)$
$\displaystyle=\left(f^{n},v_{h}\right)+\sum_{j=1}^{n-1}w_{nj}\left(-\Delta_{h}^{b_{2}}u_{h}^{j},v_{h}\right)$
(6.14)
$\displaystyle+\sum_{j=1}^{n-1}w_{nj}\left(\nabla\cdot(b_{1}(x,t_{n},t_{j})u_{h}^{j}),v_{h}\right)$
$\displaystyle+\sum_{j=1}^{n-1}w_{nj}\left(b_{0}(x,t_{n},t_{j})u_{h}^{j},v_{h}\right).$

Take $v_{h}=-\Delta_{h}^{M}u_{h}^{n}$ in (6.14) and apply (H2), (H3) along
with Cauchy-Schwarz and Young’s inequality to get

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}\nabla u_{h}^{n},\nabla
u_{h}^{n}\right)+\|\Delta_{h}^{M}u_{h}^{n}\|^{2}$
$\displaystyle\lesssim\|f^{n}\|^{2}+\sum_{j=1}^{n-1}w_{nj}\|\Delta_{h}^{b_{2}}u_{h}^{j}\|^{2}$
(6.15) $\displaystyle+\sum_{j=1}^{n-1}w_{nj}\|\nabla
u_{h}^{j}\|^{2}+\sum_{j=1}^{n-1}w_{nj}\|u_{h}^{j}\|^{2}.$

By using the estimate (4.16) and the identity
$\left(\mathbb{D}^{\alpha}_{t}u_{h}^{n},u_{h}^{n}\right)\geq\frac{1}{2}\mathbb{D}^{\alpha}_{t}\|u_{h}^{n}\|^{2}$,
the equation (6.15) is converted into to

$\mathbb{D}^{\alpha}_{t}\|\nabla
u_{h}^{n}\|^{2}+\|\Delta_{h}^{M}u_{h}^{n}\|^{2}\lesssim\left(\|f^{n}\|^{2}+\sum_{j=1}^{n-1}w_{nj}\|\nabla
u_{h}^{j}\|^{2}\right).$ (6.16)

Further proceed as we prove estimate (6.1) to complete the proof. ∎

To show the existence of the fully discrete solution $u_{h}^{n}~{}(n\geq 1)$
of the problem ($E_{\alpha}$), the following variant of Bröuwer fixed point
theorem is used.

######  Theorem 6.3.

[15] Let $H$ be finite dimensional Hilbert space. Let $G:H\rightarrow H$ be a
continuous map such that $\left(G(w),w\right)>0$ for all $w$ in $H$ with
$\|w\|=r,~{}r>0$ then there exists a $\tilde{w}$ in $H$ such that
$G(\tilde{w})=0$ and $\|\tilde{w}\|\leq r.$

######  Theorem 6.4.

Suppose that hypotheses (H1), (H2), and (H3) hold. Then there exists a unique
solution $u_{h}^{n}~{}(n\geq 1)$ to the problem ($E_{\alpha}$).

###### Proof.

(Existence) Take $n=1$ in the scheme ($E_{\alpha}$) and apply the definition
of $\mathbb{D}^{\alpha}_{t}u_{h}^{1}$ (3.4) with
$\alpha_{0}=k^{\alpha}\Gamma(2-\alpha)$ to obtain

$\displaystyle\left(u_{h}^{1}-u_{h}^{0},v_{h}\right)$
$\displaystyle+\alpha_{0}\left(M\left(x,t_{1},\|\nabla
u_{h}^{1}\|^{2}\right)\nabla u_{h}^{1},\nabla
v_{h}\right)=\alpha_{0}\left(f^{1},v_{h}\right)+\alpha_{0}kB\left(t_{1},t_{0},u_{h}^{0},v_{h}\right).$
(6.17)

In the view of (6.17) we define a map $G:X_{h}\rightarrow X_{h}$ by

$\displaystyle\left(G\left(u_{h}^{1}\right),v_{h}\right)$
$\displaystyle=\left(u_{h}^{1},v_{h}\right)-\left(u_{h}^{0},v_{h}\right)+\alpha_{0}\left(M\left(x,t_{1},\|\nabla
u_{h}^{1}\|^{2}\right)\nabla u_{h}^{1},\nabla v_{h}\right)$ (6.18)
$\displaystyle-\alpha_{0}\left(f^{1},v_{h}\right)-\alpha_{0}kB\left(t_{1},t_{0},u_{h}^{0},v_{h}\right).$

Then using (H2), (H3), and Cauchy-Schwarz inequality and Poincaré inequality
with Poincaré constant $C_{p}$, we have

$\displaystyle\left(G\left(u_{h}^{1}\right),u_{h}^{1}\right)$
$\displaystyle\geq\|u_{h}^{1}\|^{2}-\|u_{h}^{0}\|\|u_{h}^{1}\|-\alpha_{0}\|f^{1}\|\|u_{h}^{1}\|+\alpha_{0}m_{0}\|\nabla
u_{h}^{1}\|^{2}-\alpha_{0}kB_{0}\|\nabla u_{h}^{0}\|\|\nabla u_{h}^{1}\|$
(6.19)
$\displaystyle\geq\|u_{h}^{1}\|\left(\|u_{h}^{1}\|-\|u_{h}^{0}\|-\alpha_{0}\|f^{1}\|\right)+\alpha_{0}\|\nabla
u_{h}^{1}\|\left(m_{0}\|\nabla u_{h}^{1}\|-kB_{0}\|\nabla u_{h}^{0}\|\right)$
$\displaystyle\geq\|u_{h}^{1}\|\left(\|u_{h}^{1}\|-\|u_{h}^{0}\|-\alpha_{0}\|f_{h}^{1}\|\right)+\alpha_{0}\|\nabla
u_{h}^{1}\|\left(\|u_{h}^{1}\|-\frac{kB_{0}C_{p}}{m_{0}}\|\nabla
u_{h}^{0}\|\right).$

Thus, for
$\|u_{h}^{1}\|>\|u_{h}^{0}\|+\alpha_{0}\|f^{1}\|+\frac{kB_{0}C_{p}}{m_{0}}\|\nabla
u_{h}^{0}\|$ one have $\left(G\left(u_{h}^{1}\right),u_{h}^{1}\right)>0$ and
the map $G$ defined by (6.18) is continuous as a consequence of continuity of
$M$ and $B$. Hence existence of $u_{h}^{1}$ follows by Theorem 6.3
immediately.
(Uniqueness) Suppose that $X_{h}^{1}$ and $Y_{h}^{1}$ are solutions of the
scheme ($E_{\alpha}$) for $n=1$, then $Z_{h}^{1}=X_{h}^{1}-Y_{h}^{1}$
satisfies the following equation for all $v_{h}$ in $X_{h}$

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}Z_{h}^{1},v_{h}\right)$
$\displaystyle+\left(M\left(x,t_{1},\|\nabla X_{h}^{1}\|^{2}\right)\nabla
Z_{h}^{1},\nabla v_{h}\right)$ (6.20)
$\displaystyle=\left(\left[M\left(x,t_{1},\|\nabla
Y_{h}^{1}\|^{2}\right)-M\left(x,t_{1},\|\nabla
X_{h}^{1}\|^{2}\right)\right]\nabla Y_{h}^{1},\nabla v_{h}\right).$

Put $v_{h}=Z_{h}^{1}$ in (6.20) and using (H2) we get

$\frac{1}{2}\mathbb{D}^{\alpha}_{t}\|Z_{h}^{1}\|^{2}+m_{0}\|\nabla
Z_{h}^{1}\|^{2}\leq L_{M}\|\nabla Z_{h}^{1}\|\left(\|\nabla
X_{h}^{1}\|+\|\nabla Y_{h}^{1}\|\right)\left(\nabla Y_{h}^{1},\nabla
Z_{h}^{1}\right)$

Cauchy-Schwarz inequality and a priori bound (6.9) yield

$\|Z_{h}^{1}\|^{2}+k^{\alpha}\left(2m_{0}-4L_{M}K^{2}\right)\|\nabla
Z_{h}^{1}\|^{2}\leq 0.$

At last, employ (H2) to obtain $\|Z_{h}^{1}\|=\|\nabla Z_{h}^{1}\|=0$ that
concludes the uniqueness of solution for $n=1$ in the scheme ($E_{\alpha}$).

For $n\geq 2$, the numerical scheme ($E_{\alpha}$) is linear with a positive
definite coefficient matrix as a result we get the existence and uniqueness of
the solution $u_{h}^{n}~{}(n\geq 2)$ for the problem ($E_{\alpha}$). ∎

To derive the convergence rate of developed numerical scheme ($E_{\alpha}$),
first we discuss aproximation properties of L1 scheme (3.4), linearization
technique (3.5), and quadrature error (3.6).

######  Lemma 6.5.

[20] If $u\in C^{2}([0,T];L^{2}(\Omega))$ then truncation error
$\mathbb{Q}^{n}$ defined in (3.4) satisfies

$\|\mathbb{Q}^{n}\|\lesssim k^{2-\alpha}~{}\text{for}~{}n~{}\geq~{}1.$ (6.21)

######  Lemma 6.6.

For any function $u\in C^{2}[0,T]$, the linearization error
$(u^{n}-\bar{u}^{n-1}):=\bar{\mathbb{L}}^{n-1}$ defined in (3.5) undergoes

$|\bar{\mathbb{L}}^{n-1}|\lesssim k^{2}~{}\text{for}~{}n\geq 2.$ (6.22)

###### Proof.

Apply the Taylor’s series expansion of $u^{n}$ around $u^{n-1}$ and $u^{n-2}$
to obtain

$|\bar{\mathbb{L}}^{n-1}|\lesssim
k^{2}(u_{tt}(\xi_{1})+u_{tt}(\xi_{2}))~{}\text{for
some}~{}\xi_{1}~{}\in~{}(t_{n-1},t_{n})~{}\text{and for
some}~{}\xi_{2}~{}\in~{}(t_{n-2},t_{n}).$ (6.23)

As $u\in C^{2}[0,T]$ that implies the result (6.22). ∎

######  Lemma 6.7.

[29] If $u\in C^{4}[0,T]$ then quadrature error defined by (3.6) has the
following error estimate

$|q^{n}(u)|\lesssim k^{2}~{}\text{for}~{}n~{}\geq~{}1.$ (6.24)

We prove the convergence estimate of the proposed numerical scheme
($E_{\alpha}$) by assuming that the solution $u$ of the problem ($D_{\alpha}$)
satisfy the regularity assumption used in Lemma 6.5 to Lemma 6.7, i.e., $u\in
C^{4}([0,T];H^{2}(\Omega)\cap H^{1}_{0}(\Omega))$.

###  6.1 Proof of the Theorem 3.4

###### Proof.

First we prove the error estimate for the case $n=1$. Substitute
$u_{h}^{1}=W^{1}-\theta^{1}$ for $n=1$ in the scheme ($E_{\alpha}$) and using
weak formulation ($W_{\alpha}$) along with modified Ritz-Volterra projection
operator $W$ at $t_{1}$ to get

$\displaystyle\left(\mathbb{D}^{\alpha}_{t}\theta^{1},v_{h}\right)$
$\displaystyle+\left(M\left(x,t_{1},\|\nabla
u_{h}^{1}\|^{2}\right)\nabla\theta^{1},\nabla v_{h}\right)$ (6.25)
$\displaystyle=\left(\mathbb{D}^{\alpha}_{t}W^{1}-~{}^{C}D^{\alpha}_{t_{1}}u,v_{h}\right)-kB\left(t_{1},t_{0},W^{0},v_{h}\right)+\int_{0}^{t_{1}}B(t_{1},s,W(s),v_{h})ds$
$\displaystyle+\left(\left[M\left(x,t_{1},\|\nabla
u_{h}^{1}\|^{2}\right)-M\left(x,t_{1},\|\nabla u^{1}\|^{2}\right)\right]\nabla
W^{1},\nabla v_{h}\right)$
$\displaystyle+kB\left(t_{1},t_{0},\theta^{0},v_{h}\right).$

Set $v_{h}=\theta^{1}$ in (6.25) with $\theta^{0}=0$ and using (H2), (H3)
together with Cauchy-Schwarz inequality and Young’s inequality to obtain

$\displaystyle(1-k^{\alpha}\Gamma(2-\alpha))\|\theta^{1}\|^{2}+$
$\displaystyle
k^{\alpha}\left(m_{0}-4L_{M}K^{2}\right)\|\nabla\theta^{1}\|^{2}$ (6.26)
$\displaystyle\lesssim
k^{\alpha}\|\mathbb{Q}^{1}\|^{2}+k^{\alpha}\|~{}^{C}{D}^{\alpha}_{t_{1}}\rho\|^{2}+k^{\alpha}\|\nabla
q^{1}(W)\|^{2}+k^{\alpha}\|\nabla\rho^{1}\|^{2}.$

For sufficiently small $k^{\alpha}<\frac{1}{\Gamma(2-\alpha)}$, we apply (H2)
and the approximation properties (6.21), (5.5), (5.9), and (6.24) to conclude

$\|\theta^{1}\|^{2}+k^{\alpha}\|\nabla\theta^{1}\|^{2}\lesssim\left(k^{2-\alpha}+h\right)^{2}.$
(6.27)

Now we derive the error estimate for $n\geq 2$, for that take
$u_{h}^{n}=W^{n}-\theta^{n}$ in the scheme ($E_{\alpha}$)

$\displaystyle\left(\mathbb{D}^{\alpha}_{t_{n}}\theta^{n},v_{h}\right)$
$\displaystyle+\left(M\left(x,t_{n},\|\nabla\bar{u}_{h}^{n-1}\|^{2}\right)\nabla\theta^{n},\nabla
v_{h}\right)$ (6.28)
$\displaystyle=\left(\mathbb{Q}^{n},v_{h}\right)-\left({}^{C}D^{\alpha}_{t_{n}}\rho,v_{h}\right)$
$\displaystyle+\left((M(x,t_{n},\|\nabla\bar{u}_{h}^{n-1}\|^{2})-M(x,t_{n},\|\nabla
u^{n}\|^{2}))\nabla W^{n},\nabla v_{h}\right)$ $\displaystyle+\left(\nabla
q^{n}(W),\nabla
v_{h}\right)+\sum_{j=1}^{n-1}w_{nj}B(t_{n},t_{j},\theta^{j},v_{h}).$

Put $v_{h}=\theta^{n}$ in (6.28) it follows

$\displaystyle\mathbb{D}^{\alpha}_{t_{n}}\|\theta^{n}\|^{2}+\|\nabla\theta^{n}\|^{2}$
$\displaystyle\lesssim\|\mathbb{Q}^{n}\|^{2}+\|~{}^{C}D_{t_{n}}^{\alpha}\rho\|^{2}+\|\theta^{n}\|^{2}+\|\nabla\bar{\rho}^{n-1}\|^{2}+\|\nabla\bar{\theta}^{n-1}\|^{2}$
(6.29) $\displaystyle+\|\nabla\bar{\mathbb{L}}^{n-1}\|^{2}+\|\nabla
q^{n}(W)\|^{2}+\sum_{j=1}^{n-1}w_{nj}\|\nabla\theta^{j}\|^{2}.$

Employ the approximation properties (6.21), (5.5), (5.9), (6.22), and (6.24)
to deduce

$\displaystyle\mathbb{D}^{\alpha}_{t_{n}}\|\theta^{n}\|^{2}+\|\nabla\theta^{n}\|^{2}$
$\displaystyle\lesssim\left(k^{2-\alpha}+h^{2}+h+k^{2}+k^{2}\right)^{2}+\|\theta^{n}\|^{2}$
(6.30)
$\displaystyle+\left(\|\nabla\theta^{n-1}\|^{2}+\|\nabla\theta^{n-2}\|^{2}+\sum_{j=1}^{n-1}k_{1}\|\nabla\theta^{j}\|^{2}\right).$

Now follow the similar arguments as we prove estimate (6.1) to conclude the
result (3.10). ∎

##  7 Linearized L2-1σ Galerkin scheme

In this section, we show that proposed numerical scheme ($F_{\alpha}$) achieve
the second order convergence in the time direction.

######  Lemma 7.1.

Under the hypotheses (H1), (H2), and (H3) the solution $u_{h}^{n}~{}(n\geq 1)$
of the scheme ($F_{\alpha}$) satisfy the following a priori bound

$\displaystyle\max_{1\leq m\leq
N}\|u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{N}\tilde{p}_{N-n}^{(N)}\|\nabla
u_{h}^{n}\|^{2}\lesssim\max_{1\leq n\leq
N}\|f^{n-\frac{\alpha}{2}}\|^{2}+\|\nabla u_{0}\|^{2}.$ (7.1)

###### Proof.

For $n=1$ the scheme ($F_{\alpha}$) is

$\displaystyle\left(\tilde{\mathbb{D}}^{\alpha}_{t_{1-\sigma}}u_{h}^{1},v_{h}\right)+\left(M\left(x,t_{1-\sigma},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2}\right)\nabla\hat{u}_{h}^{1,\sigma},\nabla
v_{h}\right)$
$\displaystyle=\left(1-\sigma\right)kB\left(t_{1-\sigma},t_{0},u_{h}^{0},v_{h}\right)$
(7.2) $\displaystyle+\left(f^{1-\sigma},v_{h}\right).$

Substitute $v_{h}=u_{h}^{1}$ in (7.2) to get

$\displaystyle\frac{k^{-\alpha}(1-\sigma)^{1-\alpha}}{\Gamma(2-\alpha)}(u_{h}^{1}-u_{h}^{0},u_{h}^{1})$
$\displaystyle+(1-\sigma)\left(M\left(x,t_{1-\sigma},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2}\right)\nabla
u_{h}^{1},\nabla u_{h}^{1}\right)$ (7.3)
$\displaystyle=-\sigma\left(M\left(x,t_{1-\sigma},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2}\right)\nabla
u_{h}^{0},\nabla u_{h}^{1}\right)$
$\displaystyle+\left(1-\sigma\right)kB\left(t_{1-\sigma},t_{0},u_{h}^{0},u_{h}^{1}\right)+\left(f^{1-\sigma},u_{h}^{1}\right).$

Simplification of (7.3) using (H2) and (H3) with
$\alpha_{0}=k^{\alpha}\Gamma(2-\alpha)$ and
$\tilde{a}_{0}=(1-\sigma)^{1-\alpha}$ yields

$\displaystyle\left(1-\frac{\alpha_{0}}{\tilde{a}_{0}}\right)\|u_{h}^{1}\|^{2}+k^{\alpha}\|\nabla
u_{h}^{1}\|^{2}$ $\displaystyle\lesssim k^{\alpha}\|\nabla
u_{h}^{0}\|^{2}+\|u_{h}^{0}\|^{2}+k^{2+\alpha}\|\nabla
u_{h}^{0}\|^{2}+k^{\alpha}\|f^{1-\sigma}\|^{2}.$ (7.4)

Take sufficiently small $k$ to conclude

$\displaystyle\|u_{h}^{1}\|^{2}+k^{\alpha}\|\nabla
u_{h}^{1}\|^{2}\lesssim\|\nabla
u_{h}^{0}\|^{2}+\|f^{1-\sigma}\|^{2}\lesssim\|\nabla
u_{0}\|^{2}+\|f^{1-\sigma}\|^{2}.$ (7.5)

For $u_{h}^{n}~{}(n\geq 2)$ in the scheme ($F_{\alpha}$) to have

$\displaystyle\left(\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}u_{h}^{n},v_{h}\right)+\left(M\left(x,t_{n-\sigma},\|\nabla\bar{u}_{h}^{n-1,\sigma}\|^{2}\right)\nabla\hat{u}_{h}^{n,\sigma},\nabla
v_{h}\right)$
$\displaystyle=\sum_{j=1}^{n-1}\tilde{w}_{nj}B\left(t_{n-\sigma},t_{j},u_{h}^{j},v_{h}\right)$
(7.6) $\displaystyle+\left(f^{n-\sigma},v_{h}\right).$

Put $v_{h}=u_{h}^{n}$ in (7.6) to obtain

$\displaystyle\left(\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}u_{h}^{n},u_{h}^{n}\right)$
$\displaystyle+(1-\sigma)\left(M\left(x,t_{n-\sigma},\|\nabla\bar{u}_{h}^{n-1,\sigma}\|^{2}\right)\nabla
u_{h}^{n},\nabla u_{h}^{n}\right)$ (7.7)
$\displaystyle=\sum_{j=1}^{n-1}\tilde{w}_{nj}B\left(t_{n-\sigma},t_{j},u_{h}^{j},u_{h}^{n}\right)+\left(f^{n-\sigma},u_{h}^{n}\right)$
$\displaystyle-\sigma\left(M\left(x,t_{n-\sigma},\|\nabla\bar{u}_{h}^{n-1,\sigma}\|^{2}\right)\nabla
u_{h}^{n-1},\nabla u_{h}^{n}\right).$

We invoke the identity
$\left(\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}u_{h}^{n},u_{h}^{n}\right)\geq\frac{1}{2}\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}\|u_{h}^{n}\|^{2}$
and apply (H2), (H3) along with Cauchy-Schwarz and Young’s inequality to reach
at

$\displaystyle\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}\|u_{h}^{n}\|^{2}+\|\nabla
u_{h}^{n}\|^{2}$
$\displaystyle\lesssim\|f^{n-\frac{\alpha}{2}}\|^{2}+\sum_{j=1}^{n-1}\tilde{w}_{nj}\|\nabla
u_{h}^{j}\|^{2}+\|\nabla u_{h}^{n-1}\|^{2}+\|u_{h}^{n}\|^{2}.$ (7.8)

We follow the similar arguments as we prove estimate (6.1) to obtain (7.1). ∎

######  Lemma 7.2.

Under the assumptions (H1), (H2), and (H3) the solution $u_{h}^{n}~{}(n\geq
1)$ of the scheme ($F_{\alpha}$) satisfy the following a priori bound

$\displaystyle\max_{1\leq m\leq N}\|\nabla
u_{h}^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{N}\tilde{p}_{N-n}^{(N)}\|\Delta_{h}^{M}u_{h}^{n}\|^{2}\lesssim\max_{1\leq
n\leq N}\|f^{n-\frac{\alpha}{2}}\|^{2}+\|\nabla u_{0}\|^{2}.$ (7.9)

###### Proof.

We combine the idea of Lemma 6.2 and Lemma 7.1 to prove the result (7.9). ∎

######  Theorem 7.3.

Suppose that (H1), (H2), and (H3) hold. Then there exists a unique solution to
the problem ($F_{\alpha}$).

###### Proof.

For the case $n\geq 2$, the scheme ($F_{\alpha}$) is linear having positive
definite coefficient matrix. Thus existence and uniqueness in this case follow
immediately. For the case $n=1$ in the scheme ($F_{\alpha}$), the equation is
nonlinear so we again use Bröuwer fixed point Theorem 6.3. Consider the case
for $n=1$ in the problem ($F_{\alpha}$) with
$\alpha_{0}=k^{\alpha}\Gamma(2-\alpha)$ and
$\tilde{a}_{0}=(1-\sigma)^{1-\alpha}$

$\displaystyle\left(u_{h}^{1}-u_{h}^{0},v_{h}\right)$
$\displaystyle+\frac{\alpha_{0}}{\tilde{a}_{0}}\left(M\left(x,t_{1},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2}\right)\nabla\hat{u}_{h}^{1,\sigma},\nabla
v_{h}\right)$ (7.10)
$\displaystyle=\frac{\alpha_{0}}{\tilde{a}_{0}}\left(f^{1-\sigma},v_{h}\right)+\frac{\alpha_{0}}{\tilde{a}_{0}}kB\left(t_{1-\sigma},t_{0},u_{h}^{0},v_{h}\right).$

Multiply the equation (7.10) by $(1-\sigma)$ to obtain

$\displaystyle\left(\hat{u}_{h}^{1,\sigma},v_{h}\right)-\left(u_{h}^{0},v_{h}\right)$
$\displaystyle+(1-\sigma)\frac{\alpha_{0}}{\tilde{a}_{0}}\left(M\left(x,t_{1},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2}\right)\nabla\hat{u}_{h}^{1,\sigma},\nabla
v_{h}\right)$ (7.11)
$\displaystyle=(1-\sigma)\frac{\alpha_{0}}{\tilde{a}_{0}}\left(f^{1-\sigma},v_{h}\right)+(1-\sigma)\frac{\alpha_{0}}{\tilde{a}_{0}}kB\left(t_{1-\sigma},t_{0},u_{h}^{0},v_{h}\right).$

Further, proceeding analogously to the proof of Theorem 6.4 we conclude the
existence of $\hat{u}_{h}^{1,\sigma}$. Hence existence of $u_{h}^{1}$ follows.
∎

Now, we derive the convergence estimate for the numerical scheme
($F_{\alpha}$). This estimate is proved with the help of the following lemmas.

######  Lemma 7.4.

[1] If $u\in C^{3}([0,T];L^{2}(\Omega))$ then truncation error
$\tilde{\mathbb{Q}}^{n-\sigma}$ defined in (3.11) satisfies

$\|\tilde{\mathbb{Q}}^{n-\sigma}\|\lesssim
k^{3-\alpha}~{}\text{for}~{}n~{}\geq~{}1.$ (7.12)

######  Lemma 7.5.

For any function $u\in C^{2}[0,T]$, the linearization error
$(u^{n}-\bar{u}^{n-1,\sigma}):=\bar{\mathbb{L}}^{n-1,\sigma}$ defined in
(3.14) and $(u^{n}-\hat{u}^{n,\sigma}):=\hat{\mathbb{L}}^{n,\sigma}$ defined
in (3.15) undergo

$|\bar{\mathbb{L}}^{n-1,\sigma}|\lesssim k^{2}~{}\text{for}~{}n\geq 2.$ (7.13)

and

$|\hat{\mathbb{L}}^{n,\sigma}|\lesssim k^{2}~{}\text{for}~{}n\geq 1.$ (7.14)

###### Proof.

This lemma is proved by an application of Taylor’s series expansion as we have
proved Lemma (6.6). ∎

######  Lemma 7.6.

[29] If $u\in C^{4}[0,T]$ then quadrature error defined by (3.16) has the
following error estimate

$|q^{n-\sigma}(u)|\lesssim k^{2}~{}\text{for}~{}n~{}\geq~{}1.$ (7.15)

###  7.1 Proof of the Theorem 3.6

###### Proof.

Take $u_{h}^{1}=W^{1}-\theta^{1}$ with $\theta^{0}=0$ in the scheme
($F_{\alpha}$) for $n=1$, we have the following error equation for
$\theta^{1}$

$\displaystyle\left(\tilde{\mathbb{D}}^{\alpha}_{t_{1-\sigma}}\theta^{1},v_{h}\right)+\left(M\left(x,t_{1-\sigma},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2}\right)\nabla\hat{\theta}^{1,\sigma},\nabla
v_{h}\right)$ (7.16)
$\displaystyle=\left(\tilde{\mathbb{Q}}^{1-\sigma},v_{h}\right)-\left({}^{C}D^{\alpha}_{t_{1-\sigma}}\rho,v_{h}\right)+\left(M(x,t_{1-\sigma},\|\nabla
u^{1-\sigma}\|^{2})(\nabla\hat{W}^{1,\sigma}-\nabla W^{1-\sigma}),\nabla
v_{h}\right)$
$\displaystyle+\left((M(x,t_{1-\sigma},\|\nabla\hat{u}_{h}^{1,\sigma}\|^{2})-M(x,t_{1-\sigma},\|\nabla
u^{1-\sigma}\|^{2}))\nabla\hat{W}^{1,\sigma},\nabla
v_{h}\right)+\left(\nabla\tilde{q}^{1-\sigma}(W),\nabla v_{h}\right).$

Set $v_{h}=\theta^{1}$ in (7.16) with $\alpha_{0}=k^{\alpha}\Gamma(2-\alpha)$
and $\tilde{a}_{0}=(1-\sigma)^{1-\alpha}$ to have

$\displaystyle\left(1-\frac{\alpha_{0}}{\tilde{a}_{0}}\right)\|\theta^{1}\|^{2}$
$\displaystyle+k^{\alpha}\left(m_{0}-4L_{M}K^{2}\right)\|\nabla\theta^{1}\|^{2}$
(7.17) $\displaystyle\lesssim
k^{\alpha}\|\tilde{\mathbb{Q}}^{1-\sigma}\|^{2}+k^{\alpha}\|~{}^{C}D^{\alpha}_{t_{1-\sigma}}\rho\|^{2}+k^{\alpha}\|\nabla\hat{\mathbb{L}}^{1,\sigma}\|^{2}+k^{\alpha}\|\nabla\hat{\rho}^{1,\sigma}\|^{2}$
$\displaystyle+k^{\alpha}\|\nabla\tilde{q}^{1-\sigma}(W)\|^{2}.$

For small value of $k$ with hypothesis (H2) and approximation properties
(7.12), (5.5), (5.9), (7.14), and (7.15), we deduce

$\displaystyle\|\theta^{1}\|^{2}$
$\displaystyle+k^{\alpha}\|\nabla\theta^{1}\|^{2}\lesssim\left(k^{3-\alpha}+h^{2}+k^{2}+h+k^{2}\right)^{2}\lesssim\left(k^{2}+h\right)^{2}.$
(7.18)

For $n\geq 2$, substitute $u_{h}^{n}=W^{n}-\theta^{n}$ in the scheme
($F_{\alpha}$), then $\theta^{n}$ satisfies

$\displaystyle\left(\tilde{\mathbb{D}}^{\alpha}_{t_{n-\sigma}}\theta^{n},v_{h}\right)+\left(M\left(x,t_{n-\sigma},\|\nabla\bar{u}_{h}^{n-1,\sigma}\|^{2}\right)\nabla\hat{\theta}^{n,\sigma},\nabla
v_{h}\right)$ (7.19)
$\displaystyle=\left(\tilde{\mathbb{Q}}^{n-\sigma},v_{h}\right)-\left({}^{C}D^{\alpha}_{t_{n-\sigma}}\rho,v_{h}\right)+\left(M(x,t_{n-\sigma},\|\nabla\bar{u}_{h}^{n-1,\sigma}\|^{2})(\nabla\hat{W}^{n,\sigma}-\nabla
W^{n-\sigma}),\nabla v_{h}\right)$
$\displaystyle+\left((M(x,t_{n-\sigma},\|\nabla\bar{u}_{h}^{n-1,\sigma}\|^{2})-M(x,t_{n-\sigma},\|\nabla
u^{n-\sigma}\|^{2}))\nabla W^{n-\sigma},\nabla v_{h}\right)$
$\displaystyle+\left(\nabla\tilde{q}^{n-\sigma}(W),\nabla
v_{h}\right)+\sum_{j=1}^{n-1}\tilde{w}_{nj}B(t_{n-\sigma},t_{j},\theta^{j},v_{h}).$

Put $v_{h}=\theta^{n}$ in (7.19) it follows

$\displaystyle\tilde{\mathbb{D}}^{\alpha}_{t_{n-\frac{\alpha}{2}}}\|\theta^{n}\|^{2}+\|\nabla\theta^{n}\|^{2}$
$\displaystyle\lesssim\|\tilde{\mathbb{Q}}^{n-\sigma}\|^{2}+\|~{}^{C}D^{\alpha}_{t_{n-\sigma}}\rho\|^{2}+\|\nabla\hat{\mathbb{L}}^{n,\sigma}\|^{2}+\|\nabla\bar{\rho}^{n-1,\sigma}\|^{2}$
(7.20)
$\displaystyle+\|\nabla\bar{\mathbb{L}}^{n-1,\sigma}\|^{2}+\|\nabla\tilde{q}^{n-\sigma}(W)\|^{2}$
$\displaystyle+\|\nabla\theta^{n-1}\|^{2}+\|\nabla\theta^{n-2}\|^{2}+\sum_{j=1}^{n-1}\tilde{w}_{nj}\|\nabla\theta^{j}\|^{2}.$

Further, apply the approximation properties (7.12), (5.5), (5.9), (7.13),
(7.14), and (7.15) to arrive at

$\displaystyle\tilde{\mathbb{D}}^{\alpha}_{t_{n-\frac{\alpha}{2}}}\|\theta^{n}\|^{2}+\|\nabla\theta^{n}\|^{2}$
$\displaystyle\lesssim\left(\|\theta^{n}\|^{2}+\|\nabla\theta^{n-1}\|^{2}+\|\nabla\theta^{n-2}\|^{2}+\sum_{j=1}^{n-1}k_{1}\|\nabla\theta^{j}\|^{2}\right)$
(7.21)
$\displaystyle+\left(k^{3-\alpha}+h^{2}+k^{2}+h+k^{2}+k^{2}+h^{2}+h^{2}\right)^{2}.$

Now, follow the similar arguments as in Theorem 3.4 and Theorem 7.1 to obtain

$\|\theta^{m}\|^{2}+k^{\alpha}\sum_{n=1}^{m}\tilde{p}_{m-n}^{(m)}\|\nabla\theta^{n}\|^{2}\lesssim\left(h+k^{2}\right)^{2}.$

∎

##  8 Numerical results

In this section, we implement the theoretical results obtained from fully
discrete formulations ($E_{\alpha}$) and ($F_{\alpha}$) for the problem
($D_{\alpha}$). For the space discretization linear hat basis functions say
$\\{\psi_{1},\psi_{2},\dots,\psi_{J}\\}$ for $J$ dimensional subspace $X_{h}$
of $H_{0}^{1}(\Omega)$ are used, then numerical solution $u_{h}^{n}~{}(n\geq
1)$ for the considered ($D_{\alpha}$) at any time $t_{n}$ in $[0,T]$ is
written as

$u_{h}^{n}=\sum_{i=1}^{J}\alpha_{i}^{n}\psi_{i},$ (8.1)

where
$\alpha^{n}=\left(\alpha_{1}^{n},\alpha_{2}^{n},\alpha_{3}^{n},\dots,\alpha_{J}^{n}\right)$
is to be determined. Further, denote the error estimates that we have proved
in Theorems 3.4 and 3.6 by

Error-1 $\displaystyle=\max_{1\leq n\leq
N}\|u(t_{n})-u_{h}^{n}\|+\left(k^{\alpha}\sum_{n=1}^{N}p_{N-n}\|\nabla
u(t_{n})-\nabla u_{h}^{n}\|^{2}\right)^{1/2},$

and

Error-2 $\displaystyle=\max_{1\leq n\leq
N}\|u(t_{n})-u_{h}^{n}\|+\left(k^{\alpha}\sum_{n=1}^{N}\tilde{p}_{N-n}^{(N)}\|\nabla
u(t_{n})-\nabla u_{h}^{n}\|^{2}\right)^{1/2},$

respectively.

######  Example 8.1.

We consider the problem ($D_{\alpha}$) on $\Omega\times[0,T]$ where
$\Omega=[0,1]\times[0,1]$ and $T=1$ with following data

  1. 1.

$M\left(x,y,t,\|\nabla u\|^{2}\right)=a(x,y)+b(x,y)\|\nabla u\|^{2}$ with
$a(x,y)=x^{2}+y^{2}+1$ and $b(x,y)=xy$. This type of diffusion coefficient $M$
has been studied by medeiros $et.al.$ for the purpose of numerical experiments
in [24].

  2. 2.

Moreover, we take the following memory operator as in [3]
$b_{2}(t,s,x,y)=-e^{t-s}I;~{}b_{1}(t,s,x,y)=b_{0}(t,s,x,y)=0.$

  3. 3.

Source term $f(x,t)=f_{1}(x,t)+f_{2}(x,t)+f_{3}(x,t)$ with

$\displaystyle f_{1}(x,t)$
$\displaystyle=\frac{2}{\Gamma(3-\alpha)}t^{2-\alpha}(x-x^{2})(y-y^{2}),$
$\displaystyle f_{2}(x,t)$
$\displaystyle=2(x+y-x^{2}-y^{2})\left(t^{2}x^{2}+t^{2}y^{2}+\frac{xyt^{6}}{45}-2t-2+2e^{t}\right),$
$\displaystyle f_{3}(x,t)$
$\displaystyle=t^{2}(2x-1)(y-y^{2})\left(2x+\frac{yt^{4}}{45}\right)+t^{2}(2y-1)(x-x^{2})\left(2y+\frac{xt^{4}}{45}\right).$

Corresponding to the above data, the exact solution of the problem
($D_{\alpha}$) is given by $u=t^{2}\left(x-x^{2}\right)\left(y-y^{2}\right).$

We obtain errors and convergence rates in the space direction as well as in
the time direction for different parameters $h,k,$ and $\alpha$. The
convergence rate is calculated through the following $\log$ vs. $\log$ formula

$\text{ Convergence
rate}=\begin{cases}\frac{\log(E(\tau,h_{1})/E(\tau,h_{2}))}{\log(h_{1}/h_{2})};&\text{In
space direction}\\\
\frac{\log(E(\tau_{1},h)/E(\tau_{2},h))}{\log(\tau_{1}/\tau_{2})};&\text{In
time direction }\end{cases}$

where $E(\tau,h)$ denotes the error at mesh points $\tau$ and $h$.
Linearized L1 Galerkin FEM: This numerical scheme ($E_{\alpha}$) provides a
convergence order of $O(h+k^{2-\alpha})$. To observe this order of convergence
numerically, we run the MATLAB code at different iterations by setting
$h\simeq k^{2-\alpha}$. Here $h$ is taken as the area of the triangle in the
triangulation of domain $\Omega=[0,1]\times[0,1]$. For next iteration, we join
the midpoint of each edge and make another triangulation as presented in
Figures 1-3. In this way we collect the numerical results upto five iterations
to support our theoretical estimate.

Figure 1: Iteration no. 1 Figure 2: Iteration no. 2 Figure 3: Iteration no. 3

From Tables 1, 2, and 3, we conclude that the convergence rate is linear in
space and $(2-\alpha)$ in the time direction. We also observe that as
$\alpha\rightarrow 1$ then convergence rates is approaching to 1 in the time
direction which coincide with the results established in [16] for classical
diffusion case.

Table 1: _Error and Convergence Rates in space-time direction for $\alpha=0.25$ _ Iteration No. | Error-1 | Rate in Space | Rate in Time
---|---|---|---
1 | 3.04e-02 | - | -
2 | 1.50e-02 | 1.0175 | 1.7807
3 | 7.60e-03 | 0.9824 | 1.7192
4 | 3.82e-03 | 0.9909 | 1.7342
5 | 1.89e-03 | 1.0166 | 1.7790
Table 2: _Error and Convergence Rates in space-time direction for $\alpha=0.5$ _ Iteration No. | Error-1 | Rate in Space | Rate in Time
---|---|---|---
1 | 2.63e-02 | - | -
2 | 1.35e-02 | 0.9568 | 1.4352
3 | 6.53e-03 | 1.0521 | 1.5781
4 | 3.20e-03 | 1.0285 | 1.5428
5 | 1.58e-03 | 1.0141 | 1.5212
Table 3: _Error and Convergence Rates in space-time direction for $\alpha=0.75$ _ Iteration No. | Error-1 | Rate in Space | Rate in Time
---|---|---|---
1 | 2.23e-02 | - | -
2 | 1.09e-02 | 1.0340 | 1.2925
3 | 5.33e-03 | 1.0345 | 1.2931
4 | 2.65e-03 | 1.0087 | 1.2960
5 | 1.29e-03 | 1.0331 | 1.2914

Linearized L2-1σ Galerkin FEM: This numerical scheme has theoretical
convergence order of $O(h+k^{2})$. Here we set $h\simeq k^{2}$ to conclude the
convergence rates in the space-time directions. Here iteration numbers have
the same meaning as for L1 Galerkin FEM ($E_{\alpha}$). From the Tables 4, 5,
and 6, we deduce that the convergence rate is linear in space and quadratic in
the time direction.

Table 4: _Error and Convergence Rates in space-time direction for $\alpha=0.25$ _ Iteration No. | Error-2 | Rate in Space | Rate in Time
---|---|---|---
1 | 2.32e-02 | - | -
2 | 1.10e-02 | 1.0679 | 2.1359
3 | 5.38e-03 | 1.0393 | 2.0787
4 | 2.52e-03 | 1.0904 | 2.1809
5 | 1.24e-03 | 1.0225 | 2.0451
Table 5: _Error and Convergence Rates in space-time direction for $\alpha=0.5$ _ Iteration No. | Error-2 | Rate in Space | Rate in Time
---|---|---|---
1 | 2.42e-02 | - | -
2 | 1.16e-02 | 1.0573 | 2.1147
3 | 5.70e-03 | 1.0329 | 2.0658
4 | 2.68e-03 | 1.0873 | 2.1746
5 | 1.32e-03 | 1.0207 | 2.0414
Table 6: _Error and Convergence Rates in space-time direction for $\alpha=0.75$ _ Iteration No. | Error-2 | Rate in Space | Rate in Time
---|---|---|---
1 | 1.97e-02 | - | -
2 | 9.05e-03 | 1.1284 | 2.2569
3 | 4.25e-03 | 1.0873 | 2.1746
4 | 1.94e-03 | 1.1335 | 2.2671
5 | 9.30e-04 | 1.0614 | 2.1228

##  9 Conclusions

In this work, the well-posedness of the weak formulation for time-fractional
integro-differential equations of Kirchhoff type for non-homogeneous materials
is established. As a consequence of new Ritz-Volterra type projection
operator, error estimate of $O(h)$ in energy norm for semi discrete
formulation is derived. Further, to obtain the numerical solution for this
class of equations, we have developed and analyzed two different kinds of
efficient numerical schemes. First, we constructed a linearized L1 Galerkin
FEM and derived the convergence rate of order ($2$-$\alpha$) in the time
direction. Next, to enhance the convergence order in the time direction we
proposed a new linearized L2-1σ Galerkin FEM. We have proved that numerical
solutions of this scheme converge to the exact solution of the problem
($D_{\alpha}$) with the accuracy rate of $O(h+k^{2})$. Finally, numerical
results revealed that theoretical error estimates are sharp.

## References

  * [1] A. A. Alikhanov.  A new difference scheme for the time-fractional diffusion equation.  Journal of Computational Physics, 280:424–438, 2015.
  * [2] R. Almeida.  A Grönwall inequality for a general Caputo fractional operator.  Math. Inequal. Appl., 20(4):1089–1105, 2017.
  * [3] S. Barbeiro, J. A. Ferreira, and L. Pinto.  H1-second order convergent estimates for non-Fickian models.  Applied numerical mathematics, 61(2):201–215, 2011.
  * [4] J. R. Cannon and Y. Lin.  Non-classical H1 projection and Galerkin methods for nonlinear parabolic integro-differential equations.  Calcolo, 25(3):187–201, 1988.
  * [5] J. R. Cannon and Y. Lin.  A priori L2 error estimates for finite element methods for nonlinear diffusion equations with memory.  SIAM Journal on Numerical Analysis, 27(3):595–607, 1990.
  * [6] F. J. S. Correa and S. D. B. Menezes.  Existence of solutions to nonlocal and singular elliptic problems via Galerkin method.  Electronic Journal of Differential Equations, 2004:1–10, 2004.
  * [7] K. Diethelm and N. J. Ford.  Analysis of fractional differential equations.  Journal of Mathematical Analysis and Applications, 265(2):229–248, 2002.
  * [8] Y. Dimitrov.  Numerical approximations for fractional differential equations.  Journal of Fractional Calculus and Applications, 5(22):1–45, 2014\.
  * [9] Y. Giga and T. Namba.  Well-posedness of Hamilton–Jacobi equations with Caputo’s time-fractional derivative.  Communications in Partial Differential Equations, 42(7):1088–1120, 2017.
  * [10] D. Goel and K. Sreenadh.  Kirchhoff equations with Hardy–Littlewood–Sobolev critical nonlinearity.  Nonlinear Analysis, 186:162–186, 2019.
  * [11] G. Gripenberg, S. O. Londen, and O. Staffans.  Volterra integral and functional equations.  Number 34. Cambridge University Press, 1990.
  * [12] T. Gudi.  Finite element method for a nonlocal problem of Kirchhoff type.  SIAM Journal on Numerical Analysis, 50(2):657–668, 2012.
  * [13] F. Huang and F. Liu.  The time-fractional diffusion equation and the advection-dispersion equation.  The ANZIAM Journal, 46(3):317–330, 2005.
  * [14] N. I. Kavallaris, A. A. Lacey, C. V. Nikolopoulos, and C. Voong.  Behaviour of a nonlocal equation modelling linear friction welding.  IMA journal of applied mathematics, 72(5):597–616, 2007.
  * [15] S. Kesavan.  Topics in functional analysis and applications.  New Age International, New Delhi, 1989.
  * [16] L. Kumar, S. G. Sista, and K. Sreenadh.  Finite element analysis of parabolic integro-differential equations of Kirchhoff type.  Mathematical Methods in the Applied Sciences, 43(15):9129–9150, 2020.
  * [17] D. Li, H. L. Liao, W. Sun, J. Wang, and J. Zhang.  Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems.  Communication in Computational Physics, 24(1):86–103, 2018.
  * [18] L. Li and J. G. Liu.  Some compactness criteria for weak solutions of time-fractional PDEs.  SIAM Journal on Mathematical Analysis, 50(4):3963–3995, 2018.
  * [19] H. L. Liao, W. McLean, and J. Zhang.  A discrete Grönwall inequality with applications to numerical schemes for subdiffusion problems.  SIAM Journal on Numerical Analysis, 57(1):218–237, 2019.
  * [20] Y. Lin and C. Xu.  Finite difference/spectral approximations for the time-fractional diffusion equation.  Journal of Computational Physics, 225(2):1533–1552, 2007.
  * [21] J. L. Lions.  On some questions in boundary value problems of mathematical physics.  In North-Holland Mathematics Studies, volume 30, pages 284–346. Elsevier, 1978.
  * [22] F. Mainardi.  Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models.  World Scientific, 2010.
  * [23] J. Manimaran, L. Shangerganesh, A. Debbouche, and V. Antonov.  Numerical solutions for time-fractional cancer invasion system with nonlocal diffusion.  Frontiers in Physics, 7:93, 2019.
  * [24] L. A. Medeiros, T. N. Rabello, M. A. Rincon, and M. C. C. Vieira.  On perturbation of the Kirchhoff operator: analysis and numerical simulation.  Communications in Mathematical Sciences, 10(3):751–766, 2012.
  * [25] R. Metzler and J. Klafter.  The random walk’s guide to anomalous diffusion: A fractional dynamics approach.  Physics Reports, 339(1):1–77, 2000.
  * [26] K. S. Miller and B. Ross.  An introduction to the fractional calculus and fractional differential equations, John-Wily and Sons.  Inc. New York, 1993.
  * [27] X. Mingqi, V. D. Rădulescu, and B. Zhang.  A critical fractional Choquard–Kirchhoff problem with magnetic field.  Communications in Contemporary Mathematics, 21(04):1850004, 2018\.
  * [28] P. K. Mishra and K. Sreenadh.  Existence and multiplicity results for fractional p-Kirchhoff equation with sign changing nonlinearities.  Advances in Pure and Applied Mathematics, 7(2):97–114, 2016.
  * [29] A. K. Pani, V. Thomée, and L. B. Wahlbin.  Numerical methods for hyperbolic and parabolic integro-differential equations.  The Journal of Integral Equations and Applications, 4(4):533–584, 1992.
  * [30] I. Podlubny.  Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications.  Elsevier, 1998.
  * [31] W. R. Schneider and W. Wyss.  Fractional diffusion and wave equations.  Journal of Mathematical Physics, 30(1):134–144, 1989.
  * [32] V. Thomée.  Galerkin finite element methods for parabolic problems, volume 25.  Springer Science & Business Media, 2007.
  * [33] D. Tzanetis and P. M. Vlamos.  A nonlocal problem modelling ohmic heating with variable thermal conductivity.  Nonlinear analysis. Real world applications, 2(4):443–454, 2001\.
  * [34] R. Zacher.  Weak solutions of abstract evolutionary integro-differential equations in Hilbert spaces.  Funkcialaj Ekvacioj, 52(1):1–18, 2009.