A fully discrete local discontinuous Galerkin method for variable-order fourth-order equation with Caputo-Fabrizio derivative based on generalized numerical fluxes

1.   Introduction

In this paper, we will analyze the following variable-order fourth-order equations

where $0 < \alpha(t) < 1$, $\rho (t)$ is a continuous function. The solution of the problem (1.1) is periodic or compactly supported.

The variable-order Caputo-Fabrizio derivative is defined as ^[27]

Fractional calculus can better reflect the reality of nature. Many physical problems are regulated by fractional order differential equations (FDEs) and finding analytical solutions to these equations has been the subject of research by many researchers in recent years ^[5,16]. The main reason for which is that the theory of derivatives of fractions (non-integers) has aroused considerable interest in mathematics, physics, engineering and other scientific fields ^[14,21,31].

Designing an effective numerical method is meaningful for fractional order differential equations. Some numerical methods such as finite difference methods ^[4,24], orthogonal spline collocation method ^[28], finite element methods ^[18], finite volume methods ^[12], spectral methods ^[7,11], discontinuous Galerkin method ^[17,20], orthogonal spline collocation methods ^[28] and so on, have been attempted to approximate the exact solution.

Fourth-order problems as a significant part of FDEs are studied by some scholars. Combining appropriate spatial and temporal discretization, some methods have been developed to solve fourth-order FDEs, including compact difference methods ^[8,22,30], orthogonal spline collocation methods ^[23], the homotopy perturbation methods ^[3], Galerkin-Legendre spectral methods ^[2], non-polynomial quintic spline methods ^[9], finite element methods ^[13,15], LDG methods ^[6,17,25]. However, the report about numerical methods for variable-order fourth-order FDEs with Caputo-Fabrizio derivative is limited. We will study a LDG method for the problem (1.1) based on generalized numerical fluxes.

The rest of this paper is as follows. First in Section 2, some notations and necessary lemmas will be introduced. Then in Section 3, we propose a LDG scheme for solving the above problem (1.1), and discuss the stability and convergence of the method by mathematical induction. In Section 4, we shall give some numerical experiments which is made by using Matlab procedure to show the efficiency of our method. Finally in Section 5, the conclusion is given.

2.   Preliminaries

Let $a = x_{\frac{1}{2}} < x_{\frac{3}{2}} < \cdots < x_{N+\frac{1}{2}} = b$ be a partition of the domain $\bar{\Omega} = [a, b]$, $I_{j} = [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}],$ for $j = 1, \cdots N,$ and define $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\rightarrow 0^+}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)$. Futhermore, the weighted average of a function $v$ is defined by $(v)_{j+\frac{1}{2}}^{(\delta)} = {\delta}v_{j+\frac{1}{2}}^{-}+(1-\delta)v_{j+\frac{1}{2}}^{+}$, where $\delta$ is the given weight.

The local discontinuous Galerkin space $V_h^k$ is shown below

For a periodic function $\vartheta$ which is defined on the domain $[a, b]$, the generalized Gauss-Radau projections ^[1,19], denoted by ${\mathcal{P}}_{\delta}$. Let $\vartheta^e = {\mathcal{P}}_{\delta}\vartheta-\vartheta$ be the projection error. For the case $\delta \neq \frac{1}{2}$, it has the following properties

The following conclusion can be obtained from ^[1].

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

the constant $C$ which is independent of $h$, is solely dependent on the function $\vartheta$. $\tau_h$ is the union of element boundary points, and $\|\vartheta^e\|_{\tau_h}$ can be defined by

Throughout this paper, the notation $C$ represents a positive constant that may have a different value at each time. The usual notations in Sobolev space are used in the paper. Let the scalar inner product on $L^2(E)$ be denoted by $(\cdot, \cdot)_{E}$, and the associated norm by $\|\cdot\|_{E}$. If $E = \Omega$, we drop $E$.

3.   Fully discrete LDG method

Firstly, we rewrite Eq (1.1) as a system

Let $t_n = \frac{n}{M}T$, and $\tau = t_n-t_{n-1}$. The temporal derivatives $u_t$ and $\rho(t)$$_0^{CF}$$\partial_t^{1-\alpha(t)}u$ at $t_n$ are approximated as follows

where $c_k\in(t_{k-1}, t_{k})$,

and $\Phi^n(x) = \Phi^n_1(x)+\Phi^n_2(x)$ is the truncation error. According to ^[27], we can know the following conclusion

here the constant $C > 0$, depending on $T$ and the function $u$.

Furthermore, by some calculations, we could find that $b_k^n$ in Eq (3.2) has the following property

Then we could define the fully-discrete LDG method for the problem (1.1). Let $u_h^n, \; p_h^n, \; q_h^n, \; s_h^n\in V_h^k$ be the approximations of $u(\cdot, t_n), \; p(\cdot, t_n), \; q(\cdot, t_n), \; s(\cdot, t_n)$, respectively, $w^n(x) = w(x, t_n)$. Find $u_h^n, \; p_h^n, \; q_h^n, \; s_h^n\in V_h^k,$ such that for $v, \; w, \; \rho, \; \varphi\in V_h^k,$

The hat terms in Eq (3.5) which are from integration by parts are numerical fluxes. In order to guarantee stability, we will choose the following generalized numerical fluxes ^[26]

where $\delta \neq \frac{1}{2}$. If $\delta = 0$ or $1$, the flux Eq (3.6) will be purely alternating numerical fluxes ^[29].

In order to simplify the notations in the stability and convergence, we could denote

3.1.   Stability analysis

Without loss of generality, the case $w = 0$ is considered in the numerical analysis of the scheme (3.5).

Theorem 3.1. Assume that the solution of the problem (1.1) is compactly supported or periodic, then the LDG method (3.5) is stable and satisfies the following inequalities

Proof. In scheme (3.5), we first take the test functions $v = u^n_h, \; w = p^n_h, \; \rho = q^n_h, \; \varphi = -s^n_h$, we could have

In each cell $I_{j}$, we can obtain

Summing (3.10) from $1$ to $N$ over $j$, we can get

Combining Eqs (3.4), (3.11) and Cauchy-Schwarz inequality, the equality (3.9) will become

In what follows we will prove Theorem 3.1 by using mathematical induction. For the case $n = 1$ in Eq (3.12), we can easily get

that is

Next assume that the following inequality holds

we need to prove

According to Eq (3.12), we could have the following inequality

Obviously, we can directly obtain

This finishes the proof of Theorem 3.1.

3.2.   Error estimate

Theorem 3.2. Suppose that $u(x, t_n)$ is the exact solution of the problem (1.1), $u_h^n$ is the numerical solution of the fully discrete LDG scheme (3.5), then the following result holds

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

Proof.

Here $\eta_u^n$, $\eta_p^n$, $\eta_q^n$ and $\eta_s^n$ have been estimated by the inequality (2.2). Next we will estimate $\xi_u^n$, $\xi_p^n$, $\xi_q^n$ and $\xi_s^n$. Based on the fluxes (3.6), we have

Based on Eq (3.17), and taking $v = \xi_u^n, w = \xi_p^n, \rho = \xi_q^n, \varphi = -\xi_s^n$, the above error Eq (3.18) could be written as

Then, taking $v = -\xi_q^n$, $w = \xi_s^n$, $\rho = \beta\xi_u^n$, $\varphi = \beta\xi_p^n$ in Eq (3.18), we can get the following equation from the error Eq (3.18)

By applying Eqs (3.19), (3.20), stability results and $\xi_u^0 = 0$, we have the following equality

where $\; \beta = (\frac{1}{\tau}+\frac{\rho(t_n){b_n^n}}{(1-\alpha(t_n)){\tau}})$,

Next we begin to discuss the terms ${T_i}, \; i = 1, \; 2, \; 3, \; 4$.

Bsaed on the approximation property, Cauchy-Schwarz inequality, we have

and

With the help of

According to the properties (2.1), we can obtain

Noticing the fact that

and

Based on the above results, and Cauchy-Schwarz inequalities in Eq (3.21), we could have

that is

so we can obtain the following inequality

that is

Finally, by applying discrete Gronwall inequality and the triangle inequality, the proof Theorem 3.2 is completed.

4.   Numerical experiment

Consider the problem (1.1)

Let $\rho(t) = \frac{1}{2}$, and the right term $w(x, t)$ is taken such that the exact solution for the problem (1.1) is $u(x, t) = e^{t}\sin (x)$.

In the following numerical experiments, the following basis functions for $x\in\; I_j$ are choosed

The convergence results are showed for both $L^\infty$ norm and $L^2$ norm of the error. By varying the value of the parameter $\delta$, the different $\alpha(t)$, and the polynomial approximations from $P^{0}$ to $P^{2}$. For uniform meshes, numerical errors and convergence rates with different $\delta \; and\; \alpha(t)$ are shown in Tables 1–3 for $k = 0, \; 1$ and $2$, respectively. The convergent results in Tables 1 and 2 illustrate that the spatial convergence rate can attain $O(h^{k+1})$ for piecewise $P^k$ polynomials. Table 3 shows the temporal convergence rate is closed to $O(\Delta t)$ by using our scheme, this is consistent with the theoretical results.

5.   Conclusions

In this article, an accurate numerical method is presented to solve a class of variable-order (VO) fourth-order problem with Caputo-Fabrizio derivative. Based on finite difference method in time and LDG method in space, we obtain a fully discrete scheme. By taking the generalized alternating numerical fluxes, the unconditional stability and convergence are proved in detail. Some numerical experiments are shown which illustrates the effectiveness of our scheme.

Acknowledgement

This work is supported by National Natural Science Foundation of China (11861068), the Natural Science Foundation of Xinjiang Uygur Autonomous Region (2022D01E13), the Training Plan of Young Backbone Teachers in Colleges and Universities of Henan Province (2019GGJS094), Scientific and Technological Research Projects in Henan Province (212102210612), the Innovative Funds Plan of Henan University of Technology (2021ZKCJ11).

Conflicts of interest

The authors declare that they have no conflicts of interest.