A fully discrete local discontinuous Galerkin method based on generalized numerical fluxes to variable-order time-fractional reaction-diffusion problem with the Caputo fractional derivative

1.   Introduction

Fractional order partial differential equations are a generalization of classical partial differential equations ^[1,2,3]. The theory of fractional order calculus has a wide range of applications in mathematical physics equations, electrochemical processes, mechanics, anomalous diffusion, processing of signals and finance ^[4,5,6,7,8]. Since the analytical solutions of many fractional order differential equations cannot be solved exactly ^[9], numerical methods for fractional order differential equations have attracted a great deal of attention from an increasing number of scholars. Many researchers in domestic and foreign countries have solved many different types of fractional order partial differential equations by different methods. There are spectral methods ^{[10,11,12,13,14,15]}, finite difference method ^{[16,17,18,19]}, finite element method ^{[20,21,22,23,24,25]}, local discontinuous Galerkin methods (LDG) ^{[26,27,28,29,30]}. etc.

In recent years, many scholars have studied the Caputo-type reaction diffusion equation. Imran and Shah etal ^[31] investigated unsteady time fractional natural convection flow of incompressible viscous fluid, and obtained some exact solutions for temperature and velocity fields by Caputo time fractional derivatives in dimensionless form. Mahsud, Shah and Vieru ^[32] studied unsteady flows of an upper-convected Maxwell fluid which is described by the fractional differential equations with time-fractional Caputo-Fabrizio derivatives. Shah etal ^[33] studied natural convective flows of Prabhakar-like fractional viscoelastic fluids by introducing the generalized fractional constitutive equations, and used the time-fractional Prabhakar derivative to describle the generalized memory effects. J. Shu et al. ^[34] studied the asymptotic behavior of the solution of the non-autonomous fractional-order stochastic reaction-diffusion equation with multiplicative noise in R. M. Stynes et al. ^[35] studied the reaction-diffusion problem with Caputo time derivatives, giving a new analysis of the standard finite difference method for this problem. C.B. Huang et al. ^[36] presented a fully discrete numerical method for computing an approximate solution of fractional reaction-diffusion initial-boundary value problems based on $L1$ discretization in time and direct discontinuous Galerkin (DDG) finite element in space. V.K. Baranwal et al. ^[37] proposed a new analytical algorithm for solving a system of highly nonlinear time-fractional order reaction-diffusion equations, a fusion of the variational iterative method and the Adomian decomposition method. S. Ali et al. ^[38] obtained an approximate solution of the fractional order Cauchy reaction diffusion equation using the optimal homotopy asymptotic method. New numerical schemes for solving nonlinear fractional convection-diffusion equations of order $\beta\in [1,2]$ were developed by H. Safdari et al. ^[39]. They proposed locally discontinuous Galerkin methods by adopting linear, quadratic and cubic B-spline basis functions.

The Discontinuous Galerkin (DG) method is between a finite element and a finite volume method, and uses a discontinuous solution space and has high accuracy for any order of accuracy. The main idea of the LDG method is to transform the original higher-order partial differential equations into several equivalent systems of first-order equations by introducing auxiliary variables and then discretize the obtained systems of first-order equations using the DG method ^[40]. In this paper, we will consider the LDG method based on the generalized numerical flux to solve the time-fractional reaction-diffusion equation with Caputo fractional order derivatives.

in which $0 < \alpha(t)\leq \dot{\alpha} < 1$, $\sigma\geq0$. The $F, u_0$ are smooth functions. In this paper, the solution is considered to be periodic or compactly supported.

The Caputo fractional derivative operator ^[41] is defined by

In Section 2, some notations and projections are given. In Section 3, we will propose a fully discrete LDG method for the equation (1.1), and prove that the scheme is unconditionally stable and convergent with $O(h^{k+1}+(\Delta t)^{2-\dot{\alpha}})$. The correctness of the theoretical analysis is shown in Section 4 with numerical examples. Finally, the conclusion is given in Section 5.

2.   Notations and auxiliary results

Let $a = x_{\frac{1}{2}} < x_{\frac{3}{2}} < \cdots < x_{N+\frac{1}{2}} = b$ be partition of $\Omega = [a, b]$, denote $I_{j} = [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}],$ for $j = 1, \cdots N,$ and $h_j = x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}, \; 1\leq j\leq N,$ $h = \max\limits_{1\leq j\leq N}h_j$.

We denote $u_{j+\frac{1}{2}}^+ = \lim\limits_{t\rightarrow0^+}u(x_{j+\frac{1}{2}}+t)$ and $u_{j+\frac{1}{2}}^- = \lim\limits_{t\rightarrow 0^+}u(x_{j+\frac{1}{2}}-t).$

The piecewise-polynomial space $V_h^k$ is defined as

where $k$ is the order of piecewise polynomial.

For any periodic function $\varpi$, the following projection is used to prove the error estimate. Let $\omega^{e} = \mathcal{P}_{\delta}\omega-\omega$, that is $\mathcal{P}$,

Lemma 2.1. Let $\delta\neq\frac{1}{2}$. If $\omega\in H^{s+1}[a, b]$, there holds

where the bounding constant $C > 0$ is independent of $h$ and $\omega$. Here $\Gamma_h$ denotes the set of boundary points of all elements $I_j$, and

Let the scalar inner product on $L^2(E)$ be denoted by $(\cdot, \cdot)_{E}$, and the associated norm by $\Vert\cdot\Vert _{E}$. If $E = \Omega$, we drop $E$. In the paper, $C$ is a positive number that may have different values in different places.

3.   The scheme

We first construct the fully discrete LDG method for the equation (1.1).

We can rewrite the equation (1.1) into the following form:

Let $t_n = n\Delta t = \frac{n}{M}T$, $\Delta t = t_n-t_{n-1}$, we approximate the Caputo fractional derivative $_0^{C}D_t^{\alpha(t_n)}u(x, t_n)$:

where $\gamma^n$ is the truncation error in the time direction and $\Vert \gamma^n\Vert\leq C(\Delta t)^{2-\dot{\alpha}}.$

In which, $b_{k}^{n} = (k+1)^{1-\alpha(t_n)}-(k)^{1-\alpha(t_n)},$ and it can be seen that the coefficients have the following properties ^[35]

Let $u_h^n, \; p_h^n \in V_h^k$ be the approximations of $u(\cdot, t_n), \; p(\cdot, t_n)$, respectively, $F^n(x) = F(x, t_n)$. Find $u_h^n, \; p_h^n \in V_h^k$, such that for all test functions $v, w \in V_h^k,$

Selecting the appropriate numerical flux will play a key role in theoretical analysis for the LDG scheme. From the practical aspect, the generalized alternating numerical fluxes have more application than the traditional numerical fluxes ^[42]. We consider the following generalized alternating numerical fluxes

here we consider the case $\delta\neq \frac{1}{2}$. For $\delta = \frac{1}{2}$, the property about unique existence and approximation of the generalized Gauss-Radau projection will become complicated ^[43].

For the sake of convenience, we denote

Next, we give the stability analysis of the scheme (3.3).

3.1.   Stability analysis

Without loss of generality, we take $F = 0$ in the theoretical analysis. The following stability result for the scheme (3.3) can be obtained.

Theorem 3.1. For periodic or compactly supported boundary conditions, the fully-discrete LDG scheme (3.3) is unconditionally stable, and the numerical solution $u_h^n$ satisfies

Proof. Taking the test functions $v = u^n_h, \; w = p^n_h$ in the scheme (3.3), and with the fluxes choice (3.4), we obtain

In $I_{j} = [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}]$, we obtain

After some computations and summing (3.8) from $1$ to $N$ over $j$, we get

Combine (3.9) and the Cauchy-Schwarz inequality, (3.7) becomes

Use mathematical induction to prove Theorem 3.1. Let $n = 1$ in (3.10), and we can obtain

Now we assume that the following inequality holds

We need to prove

it follows from (3.10) and (3.12) that

□

3.2.   Error estimate

Theorem 3.2. Let $u(x, t_n)$ be the exact solution of the problem (1.1), which is sufficiently smooth with bounded derivatives. Let $u_h^n$ be the numerical solution of the fully discrete LDG scheme (3.3) with flux (3.4), and then the following error estimates holds

where $C$ is a constant depending on $u, T$.

Proof.

Here $\eta_u^n$ and $\eta_p^n$ have been estimated by the inequality (2.2).

Taking the flux (3.4), we can get the following error equation.

Take the test function $v = \xi_u^n, \; w = \xi_p^n$, and use (3.14) in the error equation (3.15), we can get

using the projection (2.1), and after a certain amount of analysis, we can obtain

Using the Cauchy-Schwarz inequality, we could have the following equality

let $\varpi = (\Delta t)^{\alpha(t_n)}\Gamma(2-\alpha(t_n))$, and $\varpi < \Gamma(1-\alpha(t_n))T^{\alpha(t_n)}b_{n-1}$, we have

Using the nature of projection, and noticing the fact that

we can get

so

We use mathematical induction to prove the theorem. When $n = 1$, from (3.20), it follows that

we assume that inequality holds

We only need to prove that the following results holds

Similar to the proof of stability, we can get the following formula

so

Let $\dot{\alpha} = \max\{\alpha(t_n)\}$, we have

Theorem 3.2 could be proved by using the triangle inequality and the interpolation property (2.2). □

4.   Numerical experiment

Consider the following equation (1.1)

with $u(x, 0) = 0$ for $x\in (0, 1)$, the function

is chosen such that the exact solution of the equation is $u(x, t) = t^2cos(2\pi x)$.

Choosing a fixed small time step $\Delta t = \frac{1}{1000}$ to avoid contamination of the temporal error, and dividing the space into $N$ elements to form the uniform mesh. The errors in $L^2$-norm and $L^\infty$-norm for different values $\alpha(t)$ and $\delta$s are demonstrated in Tables 1 and 2 where $N = 5, 10, 20, 40$. A uniform $(k+1)$-th order of accuracy for piecewise $P^k$ polynomials can be seen.

Finally, the temporal convergence rate of the scheme (3.3) for $\alpha(t) = \frac{1+t}{2}$ by piecewise $P^1$ polynomials are provided. Taking the spatial mesh size $h = \frac{1}{300}$, and the temporal meshes $\Delta t = \frac{1}{5}, \; \frac{1}{10}, \; \frac{1}{20}, \; \frac{1}{40}$, respectively. One can find that the temporal convergence rate is first order in Figure 1, which is also consistent with the theoretical results.

5.   Conclusions

In this paper, a numerical method is investigated to solve the variable-order (VO) fractional reaction-diffusion equation. We obtained the scheme by using the finite difference method in time and the LDG method in space. Based on the generalized alternating numerical fluxes, we prove that the method is unconditionally stable and convergent with $O(h^{k+1}+(\Delta t)^ {2-\dot{\alpha}})$. Numerical examples demonstrate the accuracy of our theoretical proofs. In the future, we will use this method to solve models for different kinds of partial differential equations.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgment

This work is supported by the Scientific and Technological Research Projects in Henan Province (212102210612).