High-order numerical method for the fractional Korteweg-de Vries equation using the discontinuous Galerkin method

1.   Introduction

Fractional differential equations have gained significant importance in recent years due to their ability to model complex phenomena such as long-range memory effects, mechanical systems, and control systems. These equations offer an enhanced framework for describing physical processes that cannot be effectively captured by integer-order models, especially in challenging scenarios like algebraic structures and noise reduction ^[1,2,3]. In particular, variable-order fractional calculus has emerged as a powerful tool in many applications, providing a natural mathematical framework for modeling systems with memory effects ^{[4,5,6,7,8,9]}.

Fractional partial differential equations (FPDEs) have demonstrated their capacity to describe anomalous physical behaviors more accurately than their integer-order counterparts, which has led to growing interest in their study. However, obtaining analytical solutions to FPDEs remains a challenge, particularly when fractional derivatives are involved. As a result, efficient numerical methods are increasingly being employed, including finite volume (element) methods ^{[10,11,12,13]}, finite difference methods ^{[15,16,17,18,19,20,21,22]}, meshless methods ^[23], spectral methods ^{[24,25,26,27]}, discontinuous Galerkin methods ^{[28,29,30,31,32,33,34,35]}, collocation methods ^[36], etc. In ^[37], Liu and Li et al. proposed a local discontinuous Galerkin (LDG) method for solving a nonlocal viscous water wave model, demonstrated stability, and provided error estimates with numerical validation. Wei and Li ^[38] proposed an exact numerical scheme for a class of variable-order fractional diffusion equations, utilizing the Caputo-Fabrizio fractional derivatives and providing a theoretical analysis using the local discontinuous Galerkin method. Du et al. ^[39] introduced various difference schemes for multi-dimensional variable-order time fractional subdiffusion equations and developed a special-point approximation for the variable-order time Caputo derivative, proving that the resulting scheme is uniquely solvable. Li et al. ^[40] conducted a numerical study on three typical Caputo-type partial differential equations using both the finite difference method and the local discontinuous Galerkin finite element method. Zhang ^[41] developed an efficient numerical scheme for solving the linearized fractional KdV equation in unbounded domains, utilizing non-local fractional derivatives, and approximating the solution of the initial-boundary value problem by exponentiating the convolution kernel.

The fractional Korteweg-de Vries (KdV) equation represents a generalization of the classical KdV equation by incorporating fractional derivatives, which allow for the modeling of truncation effects that arise in certain physical systems. Unlike the integer-order derivatives in the classical KdV equation, which primarily describe local interactions and normal dispersion, fractional derivatives capture nonlocal behaviors and memory effects. These properties are essential for modeling phenomena such as anomalous dispersion, where the wave propagation deviates from classical patterns, and long-range interactions, which are prevalent in complex systems.

In physical terms, the fractional KdV equation is particularly relevant in scenarios where the underlying dynamics involve nonlocal interactions or where the influence of past states significantly affects the present behavior. For example, in fluid dynamics, the fractional derivatives can describe wave propagation in media with complex microstructures. In plasma physics, they account for anomalous transport and long-range particle interactions. These capabilities make the fractional KdV equation a powerful tool for extending the applicability of the classical model to a broader range of real-world problems.

By incorporating these fractional effects, the fractional KdV equation provides a more accurate representation of systems where standard integer-order models fail to capture the underlying physics. In order to broaden the applicability of Korteweg-de Vries (KdV) models, it is both meaningful and challenging to develop high-order numerical schemes for solving the model equation. The discontinuous Galerkin (DG) method, which is naturally formulated for any desired order of accuracy in each element, offers flexibility and efficiency in terms of mesh construction and choice of shape functions. This method has shown great promise in various applications, including the numerical solution of nonlinear and dispersive wave equations like the KdV equation.

In this paper, we present a fully discrete LDG method based on generalized alternating numerical fluxes to solve the following fractional Korteweg-de Vries (KdV) equation:

The fractional derivative orders, denoted by $\theta(t) \in (0, 1)$, are time-dependent. The functions $f$, $\Re$, and $\varpi_0$ are assumed to be smooth, while $\varepsilon$ and $\lambda$ are positive constants. Furthermore, the solutions explored in this study are either periodic or exhibit compact support.

The Caputo-Fabrizio fractional derivative was proposed by Caputo and Fabrizio ^[42], which is defined as

The structure of this paper is as follows. In Section 2, some basic notation and mathematical foundations are introduced. Section 3 mainly introduces discrete methods and constructs the LDG scheme. Section 4 presents the stability and convergence results of the scheme. In Section 5, we give numerical experiments to illustrate the accuracy of our proposed format. Finally, we summarize and discuss our results in Section 6.

2.   Preliminaries

We divide the interval $\Omega = [a, b]$ into subintervals, denoted as $J: a = x_{\frac{1}{2}} < x_{\frac{3}{2}}< \cdots < x_{N+\frac{1}{2}} = b$. For $j = 1, \dots, N$, we define the interval $I_j = [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}]$, and the element size as $\Delta x_j = x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}$, where $h = \max_{1 \leq j \leq N} \Delta x_j$ represents the maximum element size.

We discretize the time interval $[0, T]$ uniformly into time steps, with $\Delta t = \frac{T}{M} = t_n - t_{n-1}$. The time points are given by $t_n = n \Delta t$, for $n = 0, 1, \dots, M$, which form a mesh in time.

At each boundary $x_{j + \frac{1}{2}}$, we define the left and right limits of the function $u$, denoted as $u_{j+\frac{1}{2}}^+$ and $u_{j+\frac{1}{2}}^-$, respectively. Here, $u_{j+\frac{1}{2}}^+$ refers to the value in the right cell $I_{j+1}$, and $u_{j+\frac{1}{2}}^-$ refers to the value in the left cell $I_j$. We introduce the jump notation:

The associated discontinuous Galerkin space $V_h^k$ is defined as:

where $P^k(I_j)$ represents the space of polynomials of degree $k$ over each interval $I_j$.

In this context, $C$ is a positive constant that may vary depending on the position in the expression. The notation $(\cdot, \cdot)_D$ represents the inner product over $L^2(D)$, and $\Vert \cdot \Vert_D$ represents the norm in the domain $D$. When $D = \Omega$, we omit the subscript.

3.   The LDG schemes

We bring the numerical formula for discretizing the Caputo derivative at $t = t_n$, which is expressed as ^[43]:

where

and

By further calculation we can obtain

$D_k^n$ has the following properties:

and

The truncation error in the time direction is

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

Lemma 3.1. ^[43] When $0 < \theta(t) < 1$, the truncation error ${R^n}(x)$ satisfies the following estimation:

Rewrite Eq (1.1) as a first-order system of equations:

$\varpi_h^n, \; p_h^n, \; q_h^n\in V_h^k$ represent approximate solutions of $\varpi(\cdot, t_n), \; p(\cdot, t_n), \; q(\cdot, t_n)$, respectively. $\Re^n = \Re(\cdot, t_n)$. Find $\varpi_h^n, \; p_h^n, \; q_h^n\in V_h^k$ such that for the test function $v, \; \phi, \; \chi \in V_h^k$, we have

The fluxes $\widehat{f}((\varpi_h^n)^{-}, (\varpi_h^n)^{+})$ are monotonic fluxes. Examples of monotonic fluxes suitable for local discontinuous Galerkin methods can be found ^[44]. For example, we can use the Lax-Friedrich flux, which consists of

The hat function in the element boundary term resulting from the integral by parts in (3.7) is the numerical flux. 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 ^[45]. We consider the following generalized alternating numerical fluxes:

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

4.   Stability and convergence

To streamline the notation, we focus on the scenario where $\Re = 0$ in the numerical analysis.

Theorem 4.1. Assuming periodic or compactly supported boundary conditions, the fully discrete LDG scheme (3.7) is unconditionally stable. Moreover, the numerical solution $\varpi_h^n$ satisfies the stability estimate

Proof. By selecting the test functions $v = \varpi^n_h$, $\phi = (1-\theta(t_n))\Delta t\varepsilon p_h^n$, and $\chi = -(1-\theta(t_n))\Delta t\varepsilon q_h^n$ in the discrete scheme, and applying the flux (3.8), we derive the following inequality,

Reorganizing, we obtain:

Here, the nonlinear terms are defined as:

and

Then we find:

For $n = 1$, the inequality becomes:

Using the inequality

we have:

which implies:

Suppose $\Vert \varpi_h^m \Vert \leq \Vert \varpi_h^0 \Vert$ holds for $m = 1, 2, \dots, n-1$. For $n$, from (4.3):

Dividing through by $D_n^n$, we obtain:

By induction, the unconditional stability is established. □

Next, we will derive the error estimates for the equation $f(u) = u$ under the assumption of linearity, employing (3.8) as the chosen flux formulation.

We begin by recalling some fundamental results. For any periodic function $\omega$ defined on $[a, b]$, the generalized Gauss-Radau projection ^[45], denoted by $\mathcal{P_\delta}\omega$, is the unique element in $V_h$. Let the projection error be $\omega^e = \mathcal{P_\delta}\omega - \omega$. When $\delta \neq \frac{1}{2}$, the error satisfies the following conditions for $j = 1, 2, \ldots, N$:

This result leads to the following lemma, as proven in ^[45].

Lemma 4.1. Let $\delta \neq \frac{1}{2}$. If $\omega \in H^{s+1}[a, b]$, the following inequality holds:

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

With this projection result established, we now proceed to the main error analysis.

Theorem 4.2. Let $\varpi(x, t_n)$ be the exact solution to problem (1.1), where $\varpi(x, t) \in H^{m+1}(D)$ and $0 \leq m \leq k+1$. Suppose $\varpi_h^n$ is the numerical solution obtained via the fully discrete LDG scheme (3.7). Then, the following error estimate holds:

where $C$ is a positive constant depending on $u$ and $T$.

Proof. To establish this result, we introduce the notation $\delta + \epsilon = 1$ and decompose the error as follows:

The terms $\eta_\varpi^n$, $\eta_p^n$, and $\eta_q^n$ are estimated using Lemma 4.1. Substituting these estimates into the scheme (3.7), the result follows by applying the standard stability and consistency arguments.

Based on the flux (3.8), the following error equation could be obtained:

Substituting (4.2) into (4.4), we obtain

Using the projection property (4.1) and substituting the test functions $v = \zeta_\varpi^n$, $\phi = (1-\theta(t_n)) \Delta t \varepsilon \zeta_p^n$, and $\chi = -(1-\theta(t_n)) \Delta t \varepsilon \zeta_q^n$ into (4.5), we derive the equality:

Considering that $\zeta_u^0 = 0$, we simplify the above to:

where

Applying the Cauchy–Schwarz inequality to the terms involving products, we have:

Using the inequality $ab \leq \kappa a^2 + \frac{1}{4\kappa} b^2$ with appropriate constants, we derive:

By multiplying both sides of the inequality by $2D_n^n$, we arrive at the following expression:

Using the inequality $a^2 + b^2 \leq (a + b)^2$, we obtain:

From Lemma 3.1, we know that $\| R^n \| \leq C (\Delta t)^2$ and that $(1 - \theta(t_n)) \Delta t = \mathcal{O}(\Delta t) = C \Delta t$.

Assume the following estimate holds:

When $n = 1$, we have:

which simplifies to:

Now assume that for $j = 1, 2, \dots, n-1$, the following inequality holds:

Substituting (4.12) and (4.13) into (4.10), we find:

Dividing through by $D_n^n$, we obtain:

Since $j \Delta t \leq n \Delta t = T$, we conclude:

Combining the triangle inequality and the projection property (4.1), we establish that Theorem 4.2 holds. □

5.   Numerical experiment

In this section, we analyze the effectiveness of the proposed scheme for solving the Korteweg-de Vries (KdV) equations. To illustrate the numerical performance, we consider the following example, which is defined with periodic boundary conditions and specific initial values:

where the parameters are given as $\varepsilon = 1$, $\lambda = 2$, and the nonlinear term is $f(\varpi) = \frac{1}{2} \varpi^2$. The source term is defined as

By substituting into the governing equation, it can be verified that the exact solution is given by:

This example provides a benchmark for evaluating the accuracy and stability of the numerical scheme. The exact solution serves as a reference to measure the error of the approximate solutions obtained using the proposed method.

The convergence order $s$ is defined as the rate at which the numerical error $E(h)$ decreases when the mesh size $h$ is reduced. Mathematically, the convergence order can be expressed as:

The convergence results of the proposed numerical scheme are presented in terms of both the $L^\infty$ norm and the $L^2$ norm of the error. For a uniform mesh with size $h = \frac{1}{N}$, the numerical errors and their corresponding convergence rates are summarized in Tables 1 and 2 for $k = 0$, $k = 1$, and $k = 2$, respectively. These results demonstrate the effectiveness of the method in achieving the expected convergence rates.

From the data presented in the tables, it is evident that the scheme attains the optimal convergence rates when using piecewise $P^k$ polynomials for the approximation. This aligns with the theoretical predictions and confirms the accuracy of the implementation.

From the data presented in the tables, it is evident that the scheme attains the optimal convergence rates when using piecewise $P^k$ polynomials for the approximation. This aligns with the theoretical predictions and confirms the accuracy of the implementation.

These findings highlight the flexibility and robustness of the method across different polynomial orders, confirming its capability to deliver high-order accuracy for smooth solutions. The numerical results also reinforce the theoretical claims regarding the convergence properties of the scheme.

6.   Conclusions

This paper investigates the solution of a class of time-fractional KdV equations using the local discontinuous Galerkin (LDG) method under the framework of the Caputo-Fabrizio fractional derivative. We provide a rigorous analysis of the proposed scheme, including its stability and error estimates, ensuring a solid theoretical foundation. Furthermore, numerical experiments are conducted to validate the effectiveness and demonstrate the robust numerical performance of the method.

For future work, it would be valuable to explore the extension of this scheme to more complex scenarios, such as two-dimensional or high-dimensional problems. Additionally, investigating its application to other types of fractional differential equations or real-world problems in fields like fluid dynamics, wave propagation, and material science could further establish the versatility and practicality of the proposed scheme.

Use of Generative-AI tools declaration

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

Acknowledgments

The authors would like to thank the editors and the anonymous referees for their very detailed comments and constructive suggestions, which greatly improved the presentation of this paper.

Conflict of interest

The author declares that they have no competing interests.