A new $\alpha$-robust nonlinear numerical algorithm for the time fractional nonlinear KdV equation

1.   Introduction

Korteweg-de Vries equation is a classic representative of the nonlinear dispersion equation. Since it satisfies a lot of conservation laws in solids and liquids, it is widely used in the field of gas and plasma ^[1]. Dutch mathematicians Diederik Korteweg and Gustav de Vries ^[2] first discovered the KdV equation for unidirectional motion in 1895 while studying the small and medium amplitude long wave motion of diving waves; thus, the equation is named after the above two scholars. Since then, researchers gradually found that many physical phenomena are closely related to the KdV equation such as magnetic current wave and ionic sound wave in plasma, and pressure wave in liquid-gas mixture ^[3,4].

In recent years, time fractional derivative has received extensive attention because of its heredity and memory ^{[5,6,7,8,9,10,11]}, which can simulate a large number of physical phenomena involving anomalous diffusion and non-local behaviors. There are a series of numerical works on fractional linear or nonlinear differential equations ^{[12,13,14,15,16,17,18,19,20]}. In fact, fractional derivative was first developed by pure mathematicians in the middle of the 19th century. Some 100 years later, engineers and physicists have found applications for these concepts in their areas^[21]. The emergence of the concept of fractional derivatives further strengthens the connection between disciplines. The presence of the fractional derivative introduces memory effects and non-local interactions, leading to the emergence of new wave phenomena and intriguing solutions^{[22,23,24,25]}. Correspondingly, integer-order derivatives limit the ability of equations to capture the complexity of real-world phenomena, and many phenomena involving non-local or non-local behavior cannot be accurately described using traditional integer-order derivatives. To solve these problems, one can introduce the concept of fractional derivative and apply it to the KdV equation ^[26] to obtain the following nonlinear fractional KdV equation

where $\Omega_x = (0, L), \Omega_t = (0, T]$, the source term $f(x, t)$ and initial condition $\phi(x)$ are known smooth functions. Here, the Caputo fractional derivative is used^[27,28], which is define as follows

This definition takes into account the memory effect of the system, that is, the response of the system is affected by past time periods. By introducing fractional derivatives, the nonlinear fractional derivative KdV equation can describe long-term interactions, memory effects, and non-local phenomena.

Over the past few decades, there has been an increasing amount of theoretical and numerical work on the nonlinear KdV equation. However, due to the nonlinearity and complexity of the KdV equation, it is very difficult or impossible to find the analytical solution. Therefore, a lot of scholars have devoted to obtain numerical solutions of various KdV equations. An et al.^[29] proposed a fully discrete discontinuous Galerkin (DG) method combining the well-known L1 discretization in time and DG method in space to approximate the time fractional KdV equation. Cen et al.^[30] studied a spatial first-order numerical method for integer order KdV equation with initial singularity, and a second order scheme is also presented, but no corresponding theoretical proof is given. Chen et al.^[31] studied a numerical solution of the linearized fractional order KdV equation with the initial singular on graded meshes. Shen et al. ^[32] proposed a method for the fractional KdV equation, which on graded meshes got the fact that the convergence order of the numerical scheme was $\mathcal{O}(h+N^{-\min\{2-\alpha, r\alpha\}})$. There are other KdV-types studies^{[33,34,35,36,37,38]}. According to the existing research results, there is no space second-order difference scheme with complete theory and there is much room for progress in the study of the fractional order KdV equation. Therefore, in paper we construct a spatial second-order fully discrete difference scheme with complete theory analysis. Because of the improvement of the convergence rate, it can greatly improve the operation efficiency, which is very beneficial to large-scale calculations ^[39].

In this paper, we consider an $\alpha$-robust high-order numerical method for the time fractional nonlinear KdV equation, and the major results are as follows

● We use the reduction method of to handle $u_{xxx}$ and introduce the nonlinear operator $\psi$ to approximate $uu_x$, and then a spatial second-order nonlinear difference scheme is obtained.

● We prove the existence, uniqueness, stability, and convergence of the proposed second-order nonlinear difference scheme, and then improved the unconditional stability to $\alpha$-robust stability.

● The results of numerical examples are consistent with the theory, and it is verified that the proposed nonlinear difference scheme converges with order $2$ in space and order ($2-\alpha$) in time on graded meshes.

The rest of this article is arranged as follows. In Section 2, we provide some relevant symbols and lemmas needed for theoretical proof and constructing the nonlinear difference scheme. In Section 3, we present the theoretical result, where the existence, uniqueness, unconditional stability, $\alpha$-robust stability, and convergence are proved in turn. In Section 4, two experiments are given to verify the reliability of the theoretical proof. At the end of the article, we have made a summary in Section 5.

In this paper, $C$ stands for some constants that can take on different values at different places, and the $c$ with the subscript represents a specific constant.

2.   The construction of nonlinear difference scheme

In the section, we provide some relevant symbols and lemmas needed in constructing the nonlinear difference scheme and theoretical analysis.

2.1.   Preliminaries

Given the positive integers $M$ and $N$, denote $h: = \frac L M$ be the spatial step, $x_j: = jh(0\leq j\leq M)$. Divide the interval $[0, T]$ into $N$ non-uniform compartments and set

where $r\ (r\ge 1)$ is called the grading exponent, $\tau_n = t_n-t_{n-1}$ ($n = 1, 2, \dots, N$) are time-steps. Define the grid functions as follows

Denote

and

be two set of grid functions.

Next, we introduce some important notations, for $u\in \mathcal{U}_h$, let

Assume $u, v\in \mathcal{U}_h$, we introduce inner product and norm as follows

In addition, we define the nonlinear operator $\psi$ ^[40] as follows

Now, we introduce several lemmas help to develop the theoretical analysis.

Lemma 2.1. ^[41] Let $v\in \mathcal{U}_h$ and $u\in \mathring{\mathcal{U}}_h$, then one gets

Lemma 2.2. ^[42] For any grid function $u, v\in \mathring{\mathcal{U}}_h$, then

2.2.   The construction of nonlinear difference scheme

Next, we begin to construct a nonlinear difference method for Eq.(1.1)-(1.3). Let $v = u_x$, then Eq.(1.1)–(1.3) can be written as

Now, considering Eq.(2.4a) at the points $(x_i, t_n)$ and Eq.(2.4b) at the points $(x_{i+\frac12}, t_n)$, one has

Next, we discretize the equation (2.5). First, we approximate the Caputo fractional derivative $^C_0\mathcal{D}^{\alpha}_t u(x_i, t_n)$ by employing the L1 formula on the graded meshes

where

Denote $(R_1)^n_i = ^C_0\mathcal{D}^{\alpha}_t u(x_i, t_n)-\mathit{D}^{\alpha}_N u(x_i, t_n)$, from ^[44], then we can get the error as follows

Second, using Taylor expansion and the definition of the operator $\psi$, one gets

Using second-order central difference, one arrives at

Combining Eq.(2.7)-(2.11), we can easily obtain the nonlinear difference scheme of (2.4) as follows

where the $|P^n_i| = |(R_1)^n_i+(R_2)^n_i+(R_3)^n_i|\le C(h^2+n^{-\min\{r\alpha, 2-\alpha\}}).$

Then, eliminating $P^n_i$ and $Q^n_{i+\frac12}$ in the expression and substituting numerical solution $u^n_i$ and $v^n_i$ for its exact solution $U^n_i$ and $V^n_i$, respectively, we can obtain the nonlinear difference scheme of Eq.(1.1)–(1.3) as follows

3.   Theoretical derivation

3.1.   Existence

We refer to the following Browder theorem which help us to prove the existence of the solution of the difference scheme.

Theorem 3.1. (Browder theorem)^[43] Let $\left (H, (\cdot, \cdot)\right)$ be a finite dimensional inner product space, $\|\cdot \|$ is the derived norm operator, and $\Pi:H\to H$ be continuous. Assume that

Then there exists satisfying $\left | z^*\right |\le \alpha$ such that $\Pi(z^*) = 0$.

Theorem 3.2. The nonlinear difference scheme (2.13) has at least a solution.

Proof. Denote

It is easy to get $u^0$ from (2.13c), we can get $v^0$ by computing (2.13b) and (2.13d).

Suppose $\{u^0, u^1, \cdots, u^{n-1}\}$ and $\{v^0, v^1, \cdots, v^{n-1}\}$ exist, then we now consider $\{u^n, v^n\}$ for nonlinear difference scheme (2.13), one has

and

Define $\Pi(u):\mathring{\mathcal{U}}_h \to \mathring{\mathcal{U}}_h:$

We notice that $\Pi(u)$ is a continuous operator in $\mathring{\mathcal{U}}_h$. Thereupon, taking an inner product with $u$, leads to

Then, let $\left\|u\right\| = \dfrac1{{a_{n, 1}}}(a_{n, n}\left\|u^0\right\|+\sum_{s = 1}^{n-1}(a_{n, s}-a_{n, s+1})\left\|u^{n-s}\right\|+\left\|f^n\right\|)$, we can figure out $(\Pi(u), u)\ge 0$ such that $\Pi(u^n) = 0$.

Therefore, the nonlinear difference scheme (2.13) exists a solution $\{u^n, v^n\}$ at least.

3.2.   Uniqueness

Denote:

Theorem 3.3. The solution of the difference scheme (2.13) is unique.

Proof. It is easy to get that $u^0$ and $v^0$ are unique, respectively. Now, suppose that

For $k = n$, assuming that both $\{u^n, v^n\}$ and $\{\hat{u}^n, \hat{v}^n\}$ are the solutions of (2.13), respectively, that is to say they satisfy:

and

Let

Subtracting (3.6) from (3.7), due to $u^0, u^1, \cdots, u^{n-1}$ and $v^0, v^1, \cdots, v^{n-1}$ are unique, we get

Further, taking inner product for (3.8) with $\rho^n$, one leads to

For the second term on the left side of the equation, by Lemma 2.2, leads to

Combining the previous definitions of both inner products and operators $\psi(\cdot, \cdot)$, we can derive

In addition, using Eq.(3.9) and (3.10), for the third term on the left side, we derive

Substituting (3.12)-(3.14) into (3.11), then we have

When $\frac{a_{n, 1}}{\Gamma(2-\alpha)} > \frac{c_2|\gamma|}2$, that is $\tau^{\alpha}_n < \frac{2}{\Gamma(2-\alpha)|\gamma|c_3}$. We can get $\left\|\rho^n\right\| = 0$. That is to say it holds that $u = \hat{u}, v = \hat{v}$. By mathematical induction, the solution of difference scheme (2.13) is unique.

3.3.   Stability

Theorem 3.4. ($L_2$-stability) Assume that $\{u^n_i|\ 1\le i\le M-1, 1\le n\le N\}$ is the solution of the difference scheme (2.13), then the solution is unconditionally stable.

Proof. First, we take the inner product on both sides of the Eq.(2.13a) with $u^n$, one obtains

Using Lemma 2.2, one has

Second, by (2.13b) and (2.13d), one has

Thus, we get the following inequality

Using the Cauchy-Schwarz inequality and noticing $a_{n, s}\ge a_{n, s+1}$, we arrive at

the above formula can be rearranged as

Next, we define

Then, from [Lemma 4.1 and Lemma 4.2 of ^[44]], one has

According to [Lemma 4.3, ^[44]], we select the parameter $\beta\le r\alpha$, and note that $\Gamma(1+\alpha) = \alpha\Gamma(\alpha)$, leads to

Then, we get

The above equation can also be written

Further, we can have

In the end, combined Eq.(3.23) and Eq.(3.27), we get

Remark 3.5. In this theorem, we have completed the stability proof of solution. However, when $\alpha\to1^-$, we note that it leads to $\Gamma(1-\alpha)\to \infty$, and the values on the right-hand side of the inequality are no longer binding. In order to avoid this shortcoming, we next to improved the result in Theorem 3.7.

Lemma 3.6. ^[45] For any finite time $t_N = T > 0$ and a given nonnegative sequence $(\lambda_l)^{N-1}_{l = 0}$, assume that there exists a constant $\lambda$, independent of time-steps, such that $\lambda\ge \sum_{l = 0}^{N-1}\lambda_l$. Suppose that the grid function $\{v^n|n\ge 0\}$ satisfies

where $\{\xi^n, \eta^n|1\le n\le N\}$ are nonnegative sequences. If the time-step satisfies $\tau_{k-1}\le \tau_k$ and the maximum time-step $\tau_N\le \sqrt[\alpha]{\frac{1}{2\Gamma(2-\alpha)\lambda}}$, one holds that

where $E_{\alpha}(x) = \sum_{k = 0}^{\infty}x^k/\Gamma(1+k\alpha)$ is the Mittag-Leffler function, $\omega_{\beta}(t) = \frac{t^{\beta-1}}{\Gamma(\beta)}$, and the discrete convolution kernel $P^{(n)}_{n-k}$ is defined as follows

Theorem 3.7. ($\alpha$-robust $L_2$-stability) Assume that $\{u^n_i|1\le i\le M-1, 1\le n\le N\}$ is the solution of the scheme (2.13), then one has

Proof. First, taking inner product to $u^n$ for the Eq.(2.13a), we can get

According to Lemma 2.2 and equation (3.18), we have

and

Therefore, employing Cauchy-Schwartz inequality and Young inequality, we arrive at

Further, due to the Lemma 3.6, one leads to

Finally, we can get the result of stability

which is equivalent to

The above theorem solves the shortcoming when $\alpha\to1^-$.

3.4.   Convergence

In this subsection, we will prove the convergence of the nonlinear difference scheme (2.13). First, we introduce the following lemma to help us complete the proof.

Lemma 3.8. ^[46] For any fixed $t\in [0, T]$, if $u(x, \cdot)\in C^6([0, L])$, for $0\le k \le N$, denote

Then there exists a constant $c_3$, one has

and

Theorem 3.9. Assume $\{U^k_i, V^k_i|\ 0\le i\le M, \ 0\le k\le N\}$ and $\{u^k_i, v^k_i|\ 0\le i\le M, \ 0\le k\le N\}$ be the solution of (2.12) and the difference scheme (2.13), respectively. Then there exists a constant such that

Proof. Denote

and

Subtracting (2.12) from (2.13), one can get the system of error equation

Taking inner product of (3.35) with $e^n$, leads to

For the first term, by the definition and Cauchy-Schwartz inequality, we have

According to equation (3.13), for the second term on right of equation, one has

For the second term on the left hand

Rewrite $g^n_i$ as follows

By the definitions of $R^n_i$ and $\delta_xR^n_{i+\frac12}$, we arrive at

Substituting the above results and applying Lemma 3.8, leads to

Combining with (3.39)-(3.45), we get

Further, there exists a constant $c_5$ such that

Using the Lemma 3.6, then we have

Taking the square root of both sides, for simplicity, the above equation can be written as

Multiplying both sides of this inequality by $\tau_n$, and then summing up for $n$ from $1$ to $N$, arrives at

We select $r = \frac{2-\alpha}{\alpha}$, then we can get the optimal grids, which leads to the difference scheme can achieve ($2-\alpha$) precision in the time direction.

4.   Numerical Experiments

In this section, we will present two numerical experiments to verify the reliability of the previous theoretical results. From the theoretical proofs of the nonlinear difference scheme established in this paper, it can be found that the scheme has only one solution, so we choose the fixed point iteration method to calculate the numerical solution. For details, see reference ^[47]HY__HY, Algorithm 1]. Another way, for the variable $v$ introduced in the nonlinear difference scheme, it is not necessary to participate in the calculation, and the difference scheme can be separated to obtain a difference system with only variable $u$. We select $T = L = \gamma = 1$ and employ the scheme (2.13) to solve the following examples. In order to achieve the optimal convergence rate in time, we select $r = \frac{2-\alpha}{\alpha}$.

In our examples we consider the exact solution has no known and define $L_{\infty}$-errors and convergence rate, respectively.

where the $u$ is reference solution that can approximately replace the exact solution in time or space with high-precision.

Example 1. In this example, we consider the initial value condition is $u^0(x) = 0$ and the source term is

In Example 1, we can compute numerical solutions for different nodes at different time and space steps as in Table 1. To verify the spatial convergence rate, we first fix $N = 512$. Table 2 lists the convergence rates for different $\alpha$. We can observe that the spatial convergence rate is order $2$, and it is clear that the calculated result satisfies the theoretical expectation well. On the other hand, when we verify the temporal convergence rate, we fix $M = 512$. Table 3 lists the change of the numerical solution of Example 1 with $N$ for different $\alpha$, where $Rate_t$ is the convergence rate in the time direction. However, we notice that both $\alpha$ and $N$ values are small, the convergence rate is slow, failing to meet the theoretical expectation of ($2-\alpha$). Thus, we increase the number of nodes, and the convergence rate gradually approaches the theoretical value. Therefore, the time convergence diagram drawn in Figure 3, the space convergence diagram drawn in Figure 4, where the numerical curve are parallel to the theoretical expected curve. This shows the superiority of graded meshes.

Example 2. For the second example, we consider the initial condition is $u^0(x) = 0$, and the source term is

In Example 2, fixed $N = 512$, Table 2 lists the changes of the maximum error values $E(M, N)$ of different $\alpha$ with $M$. From the table, we see that the spatial convergence rate is order $2$. Then, to verify the time convergence rate, fixed $M = 512$, Table 5 lists the variation of the maximum error $E(M, N)$ with $N$ for different $\alpha$. Similar to Example 1, when $\alpha$ is larger, the time convergence rate is ($2-\alpha$), which is consistent with the theory. When $\alpha$ is small, due to the effect of the singularity of this equation, the time convergence rate is lower than the theoretical value, so we increase the time node quantity $N$, and with the increase of $N$, the time convergence rate gradually returns to ($2-\alpha$).

The solution of the KdV equation describes the propagation and evolution behavior of waves in the system under the influence of nonlinear, non-local and memory effects. The dynamic evolution graph enables us to understand the spatial and temporal characteristics of the solution, the physical meaning of the evolutionary behavior, and the propagation and interaction of the wave. Figures 1–2 and Figures 5–6 are the dynamic evolution graph of the solutions of Example 1 and Example 2 respectively. Because the temporal convergence rate is ($2-\alpha$), the numerical solution obtained by taking different $\alpha$ at the same time will change with the change of $\alpha$. Figures 1 and Figure 2 draw the numerical solution surface obtained by different $\alpha$ of Example 1, Figure 5 and Figure 6 draw the numerical solution surface obtained by different $\alpha$ of Example 2, from which we can see only the slight difference. To highlight the difference, we draw Figure 7, the curve of the value obtained by different $\alpha$ at the same time $T = 1$ for the example 2.

5.   Summary

In this paper, a nonlinear difference scheme has been established for solving time fractional nonlinear KdV equations. The precision of $\tau^{2-\alpha}$ is achieved in time and $h^2$ in space. In theoretical analysis, the mathematical induction is applied in the proof of the existence and uniqueness, and the Browder theorem plays an important role for the existence analysis of nonlinear difference scheme solutions. In addition, the unconditional stability proof of the nonlinear difference scheme and the $\alpha$-robust stability proof are given. Finally, the proof of convergence in the $L_2$ norm is derived. The above several theoretical results are verified by two numerical examples. It is worth noting that two kinds of stability analysis are carried out and that the first stability result will explode when $\alpha \to 1^-$, so we introduce the discrete Gronwall inequality to improve the analysis. The relevant theoretical proofs in this paper are complementary to the fractional order nonlinear KdV equation.

Use of AI tools declaration

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

Acknowledgments

We thank the anonymous referees for their valuable comments and suggestions which helped us to improve the manuscript a lot.

The work was supported by National Natural Science Foundation of China Mathematics Tianyuan Foundation (12226337, 12226340, 12126321, 12126307), Scientific Research Fund of Hunan Provincial Education Department (21B0550, 22C0323, 23C0193), Hunan Provincial Natural Science Foundation of China (2022JJ50083, 2023JJ50164).

Conflict of interest

The authors declare that there no conflicts of interest.