Numerical analysis of variable-order fractional KdV-Burgers-Kuramoto equation

1.   Introduction

Fractional calculus is the development and extension of integral calculus, and its practical significance has attracted wide attention of many scholars ^[1,2,3]. Its theoretical research has been well applied not only in the field of pure mathematical theory, but also in various fields such as rheology and mechanical systems, optical and thermal systems, signal and image processing, industrial production and technology, electromagnetism, physics and materials ^{[4,5,6,7,8,9]}. In recent years, variable order fractional calculus has been discovered in some physical processes, such as algebraic structure and noise reduction. It is a generalization of the concepts of constant order fractional differential and fractional differential theory ^{[10,11,12,13,14,15,16,17,18]}.

KdV-Burgers-Kuramoto equation is a partial differential equation in nonlinear mathematical physics equation. The equation includes unsteady term, nonlinear term, dissipation term, dispersion term and instability term ^[19]. KdV class equations are the classical equations of nonlinear wave theory and soliton phenomenon. It is first obtained from the shallow water wave equation, and appears in plasma acoustic waves, longitudinal dispersion waves in elastic rods, etc., which can well explain many important physical phenomena. It can be used to describe the movement of unstable systems such as turbulence in some physical processes, unstable floating waves in plasma, stress waves in broken porous media, and nonlinear surface long waves of viscous fluids flowing downward on inclined surfaces ^[20,21]. Some papers have studied nonlinear or linear differential equations using local discontinuous Galerkin or other methods ^{[22,23,24,25,26,27]}, while others have studied fourth-order differential equations ^{[28,29,30,31,32]}. In recent years, many scholars have studied various applications and properties of KdV-Burgers-Kuramoto equation and different numerical methods. Secer and Ozdemir proposed the wavelet Galerkin method to solve the time fractional order KdV-Burgers-Kuramoto equation, which was transformed into its corresponding nonlinear algebraic equation, and the numerical solution was obtained by Newton method ^[33]. Bruzón et al. introduced the generalized KdV-Burgers-Kuramoto equation, and studied the conservation law and classical symmetry with multiplier method ^[34]. Kim and Chun found the exact solution of the equation based on empirical function method, and constructed a new generalized solitary wave solution of KdV-Burgers-Kuramoto equation by using truncated Painleve expansion method and Exp function method ^[35]. Kaya et al. linearized the nonlinear term of the equation by linearization technique, and obtained its numerical solution by finite difference method ^[36].

In this paper, we will design a high-order local discontinuous Galerkin method to simulate the fractional KdV-Burgers-Kuramoto equation with variable-order time derivatives, and discuss its stability and the optimal convergence rate. Discontinuous Galerkin (DG) method not only uses the element polynomial space of general finite element method as the approximate solution and test function space, but also allows the basis function to adopt completely discontinuous piecewise polynomial, and has the advantage of both finite element method and finite volume method ^[37,38]. DG plays an important role in volume dynamics, shallow water simulation, magnetohydrodynamics, oceanography, viscoelastic flow, oil recovery simulation, semiconductor device simulation and so on ^{[39,40,41,42]}.

We consider the following variable-order time fractional KdV-Burgers-Kuramoto equation

in which $0 < \epsilon(t) < 1$, $\theta_1$, $\theta_2$, $\theta_3$ are constants and $\theta_1, \theta_3\geq0$. The $L(u)$ is an arbitrary nonlinear function and $F, u_0$ are smooth functions. In this paper, the solution is considered to be periodic or compactly supported.

The variable-order Riemann-Liouville fractional derivative in Eq (1.1) is defined by

In Section 2, some symbols, basic projections and the numerical flux are given. In Section 3, we will propose a fully discrete local discontinuous Galerkin method for the Eq (1.1), and prove that the scheme is unconditional stable and convergent. Numerical examples are given to show the reliability and effectiveness of the method in Section 4. Finally, the conclusion is given in Section 5.

2.   Notations and auxiliary results

2.1.   Notations and projection

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).$

Let $[u_h^n]_{j+\frac{1}{2}}$ is used to denote $(u_h^n)_{j+\frac{1}{2}}^+- (u_h^n)_{j+\frac{1}{2}}^-,$ the jump of $u_h^n$ at each element boundary point.

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

where $k$ is order of piecewise polynomial.

The fractional derivative of Riemann-Liouville and Caputo are introduced below, which are related to each other ^[43,44].

Lemma 2.1. Let $W^m_p(a, b)$ be the Banach spaces with their weak derivatives of order $m$ in $L^p(a, b)$, and the function $\mu(x, t)\in W^1_1(0, T)$, then we have

where $_0^{C}D_t^{1-\epsilon(t)}\mu(x, t)$ is the variable-order Caputo fractional derivative

For any periodic function $\varpi$, the following are two basic projections that will be used in error analysis, that is $\mathcal{P}$,

and projection $\mathcal{P}^\pm$,

and

The projections $\mathcal{P}$ and $\mathcal{P}^\pm$ satisfy^[45,46,47]

where $v = \mathcal{P}\varpi-\varpi$ or $v = \mathcal{P}^\pm\varpi-\varpi$. The positive constant $C$, solely depending on $\varpi$, is independent of $h$. $\tau_h$ represents the union of all element boundary points, and the $L^2$-norm on $\tau_h$ is defined by

2.2.   Numerical flux

Numerical flux $\widehat{L}(\psi^-, \psi^+)$ is considered in this paper. It is monotone, which depends on the two values of the function $\psi$ at the discontinuity point $x_{j+\frac{1}{2}}$. Many examples of monotonic flux can be found in reference ^[48].

In the paper, $C$ is a positive number that may have different values in different places. 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$.

3.   Fully discrete LDG scheme

We first describe the fully discrete local discontinuous Galerkin method for the Eq (1.1). By means of Lemma 2.1, we can rewrite the model Eq (1.1) into the following form

Let $t_n = n\Delta t = \frac{n}{M}T$, $\Delta t = t_n-t_{n-1}$. We estimate the time derivative $_0^{C}D_t^{1-\epsilon(t)}u$ at $t_n$ as follows

where $\omega_i^n = (i+1)^{\epsilon(t_n)}-i^{\epsilon(t_n)}$. The truncation error is

Further, we have the following results:

where $\omega_i^n$ has the following properties

Let $u_h^n, \; g_h^n, \; \iota_h^n, \; \gamma_h^n, \; \beta_h^n \in V_h^k$ be the approximations of $u(\cdot, t_n), \; g(\cdot, t_n), \; \iota(\cdot, t_n), \; \gamma(\cdot, t_n), \; \beta(\cdot, t_n)$, respectively, $F^n(x) = F(x, t_n)$. Find $u_h^n, \; g_h^n, \; \iota_h^n, \; \gamma_h^n, \; \beta_h^n \in V_h^k$, such that for all test functions $\varrho, \; \sigma, \; \varsigma, \; \kappa, \; \chi \in V_h^k,$

For the sake of convenience, we take $\lambda = \frac{\Gamma({1+\epsilon(t_n)})}{\Delta t^{\epsilon(t_n)-1}}$. The choice of the hat term (3.6) will have an important influence on the stability. We consider the following numerical flux

where $\tau > 0$.

We use iterative method to calculate it because of the nonlinear.

Find $u_h^{n, m}, \; g_h^{n, m}, \; \iota_h^{n, m}, \; \gamma_h^{n, m}, \; \beta_h^{n, m} \in V_h^k$, such that for all test functions $\varrho, \; \sigma, \; \varsigma, \; \kappa, \; \chi \in V_h^k$, we can get

Where $m$ is the iteration step. $u_h^{n, 0} = u_h^{n-1}$ is initial condition, $\Vert u_h^{n, m}-u_h^{n, m-1}\Vert\leq10^{-6}$ is stop condition.

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

3.1.   Stability analysis

Without losing generality, we consider the case of $F = 0$ and (3.8) in the numerical analysis of this model problem. The following stability result for the scheme (3.6) is obtained.

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

Proof. Taking the test functions $\varrho = u^n_h, \; \chi = -\lambda g^n_h+\lambda\theta_1 \beta^n_h, \; \sigma = \lambda \beta^n_h, \; \chi = \lambda(\iota^n_h+\theta_1 u^n_h+|\theta_2|\beta^n_h), \; \varsigma = -\lambda \gamma^n_h$ in scheme (3.6), and with the fluxes choice (3.7), we obtain

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

After some calculation, and sum (3.12) from $1$ to $N$ over $j$, we can easily get

Let $\dot{L}(u) = \int_{0}^{u}L(u)du$. Considering the nonlinear term, we can use a mean value theorem to calculate, so that we can get

where $\xi$ is a value between $(u_h^n)^{-}$ and $(u_h^n)^{+}$. We consider the monotonicity of flux function, and obtain inequality. Combine (3.5) and (3.13), the equality (3.11) becomes

The Theorem 3.1 will be proved by mathematical induction. Let $n = 1$ in (3.14), and based on the following formula:

We can obtain

then we can get the following inequalities immediately

and

Now we assume that the following inequality holds

we need to prove

It follows from (3.14) that

Consequently, we have

3.2.   Error estimate

Consider the linear case $L(u) = u$ and choose (3.8) as the numerical flux in error estimate. We have the following theorem.

Theorem 3.2. Let $u(x, t_n)$ be the exact solution of theproblem (1.1), which is sufficiently smooth with boundedderivatives. Let $u_h^n$ be the numerical solution of the fullydiscrete LDG scheme (3.6), then there holds the followingerror estimates when $0 < \epsilon(t_n)\leq\bar{\epsilon} < 1$

and when $\epsilon(t_n)\leq\bar{\epsilon}\rightarrow 1$,

Proof.

Here $\eta_u^n$, $\eta_g^n$, $\eta_\iota^n$, $\eta_\gamma^n$, and $\eta_\beta^n$ have been estimated by the inequality (Eq 2.4).

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

Take the test function $\varrho = \xi_u^n, \; \chi = -\lambda\xi_g^n+\lambda\theta_1\xi_\beta^n, \; \sigma = \lambda\xi_\beta^n, \; \kappa = \lambda(\xi_\iota^n+\theta_1\xi_u^n+|\theta_2|\xi_\beta^n), \; \varsigma = -\lambda\xi_\gamma^n$, and use (3.20) in the error equation (3.21), we can get

From the stability result (3.13), and notice that $\omega_{n-1}^n > \epsilon(t_n)n^{\epsilon(t_n)-1}$, we could have the following equality

Noticing the fact that $ab\leq\varepsilon a^2+\frac{1}{4\varepsilon}b^2$, using the Hold's inequality, we have

choosing a small enough $\varepsilon$, we have

The error estimation will be proved by mathematical induction. For the sake of convenience, we denote

1) We assume that inequality holds

When $n = 1$, the Eq (3.25) becomes

it is easy to see that $\Vert \xi_u^0\Vert \leq C_0h^{k+1}$, we use the projection (2.4), the Eq (3.27) becomes

denoting $C = \max\{3C_0, C_1C_2, \sqrt{(\frac{C_2(1+(\tau|\theta_2|)^2+2\varepsilon+\theta_1^2\varepsilon)}{\varepsilon})}C_0 \}$, then we can obtain

When $n = k+1$, from the Eq (3.25), we can get the following formula

from the fact that $\omega_{i-1}^{-1} < \omega_{i}^{-1}$, we can obtain

that is

The inequality (Eq 3.26) follows.

By some calculations and analysis, we know that $n^{-\epsilon(t_n)}\omega_{n-1}^{-1}$ increasingly tends to $\frac{1}{1-\epsilon(t_n)}$. So we can obtain

let $\bar{\epsilon} = \max\{\epsilon(t_n)\}$, we get

2) The above estimate has no any meaning when $\epsilon(t_n)\rightarrow 1$ due to $\frac{1}{1-\epsilon(t_n)}\rightarrow \infty$. So we must reconsider it for the case $\epsilon(t_n)\rightarrow 1$. We suppose the following estimate holds

By the similar techniques used in 1), we can obtain (3.33) easily. Here we omitted the proof to save space. Then we know that when $\epsilon(t_n)\rightarrow 1$,

Therefore, the Theorem (3.2) is proved by using triangle inequality and interpolation property (2.4).

4.   Numerical experiment

Consider the following Eq (1.1)

with $u(x, 0) = 0$ for $x\in (0, 2\pi)$. Let $L(u) = \frac{1}{2}u^2$, and the function

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

The convergence results are obtained for both $L^\infty$ norm and $L^2$ norm of the error. For uniform meshes of size $h = \frac{1}{N}$, numerical errors and convergence rates are shown in Table 1 for $k = 0, \; 1$ and $2$, respectively. The approximate results illustrate that we can obtain the optimal convergence rate for piecewise $P^k$ polynomials.

5.   Conclusions

In this paper, a fully discrete local discontinuous Galerkin finite element method for solving the nonlinear variable order KdV-Burgers-Kuramoto equation is presented, which based on the finite difference method and the local discontinuous Galerkin method. By choosing the numerical flux carefully, we prove that the scheme is unconditionally stable and convergent. Numerical results show that the method is effective for solving this kind of equations. In the future, we will develop the method discussed in this paper to solve various fractional problems in physical processes.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (12126315, 12126325, 12026263, 12026254), the Fundamental Research Funds for the Henan Provincial Colleges and Universities in Henan University of Technology, PR China (2018RCJH10), the Training Plan of Young Backbone Teachers in Henan University of Technology, PR China (21420049), the Training Plan of Young Backbone Teachers in Colleges and Universities of Henan Province, PR China (2019GGJS094), Scientific and Technological Research Projects in Henan Province (212102210612), and the Humanities and Social Science Youth Foundation of Ministry of Education of China (17YJC630084).

Conflict of interest

The authors declare there is no conflicts of interest.