Non-Standard Finite Difference and Vieta-Lucas Orthogonal Polynomials for the Multi-Space Fractional-Order Coupled Korteweg-de Vries Equation

Abstract

This paper focuses on examining numerical solutions for fractional-order models within the context of the coupled multi-space Korteweg-de Vries problem (CMSKDV). Different types of kernels, including Liouville-Caputo fractional derivative, as well as Caputo-Fabrizio and Atangana-Baleanu fractional derivatives, are utilized in the examination. For this purpose, the nonstandard finite difference method and spectral collocation method with the properties of the Shifted Vieta-Lucas orthogonal polynomials are employed for converting these models into a system of algebraic equations. The Newton-Raphson technique is then applied to solve these algebraic equations. Since there is no exact solution for non-integer order, we use the absolute two-step error to verify the accuracy of the proposed numerical results.

1. Introduction

Mathematicians who delve into fractional calculus broaden the principles of differentiation and integration to include orders that are non integer. Fractional calculus allows for differentiation and integration to be performed using non-integer, rather than being confined to whole numbers or integers [1,2,3,4]. In the 17th century, pioneers like Isaac Newton and Gottfried Leibniz made early strides in advancing fractional calculus. However, it was not until the 19th century that mathematicians such as Karl Weierstrass and Augustin-Louis Cauchy embarked on the development of a structured theory for fractional calculus [3].

There are alternative approaches to defining fractional derivatives and integrals, which encompass the Riemann-Liouville, Liouville-Caputo, and Grünwald-Letnikov definitions. These definitions entail expanding the classical methods of integration and differentiation to include orders that are not integer [5].

Fractional calculus has extensive use across various industries like engineering, physics, finance, and signal processing. It offers an advanced approach to comprehend phenomena like as anomalous diffusion, fractal time series, and, viscoelasticity, all of which exhibit fractional or fractal characteristics [6,7,8,9,10,11,12,13,14]. In addition, fractional calculus provides a potent mathematical tool for modeling and examining complex systems that cannot be entirely elucidated by employing the standard methods of integer-order calculus [15].

When studying nonlinear dispersive waves, researchers often turn to the Korteweg-de Vries (KdV) equations. For more information [16].

In 1895, Korteweg and de Vries formulated these equations with the aim of simulating shallow water waves in a canal [17].

A nonlinear partial differential equation known as the coupled Korteweg-de Vries (cKdV) equation describes the evolution of many broad, weakly nonlinear waves in a dispersed medium. It is a generalization of the well-known Korteweg-de Vries (KdV) equation, which describes the behavior of an isolated single wave in particular. Hirota and Satsuma first introduced the development of the linked KdV equation in a study identified as [18].

The KdV equations elucidate the interaction between two lengthy waves characterized by distinct dispersion relations [18,19].

It is important to highlight that the Korteweg-de Vries equations appear in diverse fields within the physical sciences, ranging from plasmas and fluids to vibrations in low-temperature crystal lattices. Although these applications may appear unrelated, they all originate from a comprehensive physical model. When a particular aspect of the problem is examined closely, they all lead to the Korteweg-de Vries equation. In this sense, the Korteweg-de Vries equation is considered to have broad and universal applicability.

Vieta-Lucas polynomials are widely used in spectral methods used to solve partial and ordinary differential equations [20]. These methods show fast exponential convergence or spectral precision, particularly for problems in simple geometric contexts where the solutions are smooth. These methods differ from finite difference techniques in that the determination of the approximation coefficient entirely defines the solution at any specified point within the interval. Therefore, utilizing operational matrix spectral methods that rely on Vieta-Lucas polynomials for solving the system outlined in Equation (1) below is of considerable significance [21,22,23]. For the application of spectral collocation method and finite difference method, some powerful mathematical techniques have been used in recent years by researchers [24,25,26,27,28].

The objective of this study is to offer numerical solutions for the multi-space fractional-order coupled Korteweg-de Vries problem. Various kernels, Liouville-Caputo fractional derivatives, Caputo-Fabrizio, and Atangana-Baleanu fractional derivatives are employed in this work. The solution method integrates the Vieta-Lucas collocation method with the nonstandard finite difference method to discretize the spatial fractional-order diffusion equation. This results in a set of ordinary differential equations, simplifying the conditions significantly. The nonstandard finite difference method (FDM) is then applied to solve this set of ordinary differential equations. For more detailed information, refer to the provided references [29,30,31].

The paper is organized as follows: In Section 2, we delve into function approximation and the application of Vieta-Lucas polynomials. Section 3 demonstrates the use of the spectral method to solve the multi-space fractional-order coupled Korteweg-de Vries problem, employing Liouville-Caputo, Caputo-Fabrizio, and Atangana-Baleanu fractional derivatives. Section 4 presents an explanation and numerical results for the multi-space fractional-order coupled Korteweg-de Vries, incorporating the three kernels utilized in this study. Finally, concluding remarks are provided in Section 5.

2. Vieta-Lucas Polynomials and Function Approximations

In this section, we furnish the definitions of the shifted Vieta-Lucas polynomials (VLPs), along with their corresponding notations and key properties. This serves as a foundational understanding for their application in the subsequent analyses [32]. Our research has primarily focused on a specific category of orthogonal polynomials known as the Vieta-Lucas polynomials. This family of orthogonal polynomials can be generated through the utilization of recurrence relations and analytical equations associated with these polynomials. Also we will introduce the convergence analysis and the nonstandard finite difference scheme notations.

2.1. Shifted Vieta-Lucas Polynomials

The Vieta-Lucas polynomials ${\mathsf{\Psi }}_j\left(z\right)$, referred to as Lucas polynomials (VLPs), are discussed along with their respective notations and properties as outlined by Vieta. Specifically, for polynomials of degree $j\in {\mathbb{N}}_0$, ${\mathbb{N}}_0=\{0,1,2,3\dots \},$ the definition is provided as follows [32]:

$${\mathsf{\Psi }}_j\left(t\right)=2cos\left(j\psi \right),\psi =arccos\left(0.5z\right),\psi \in [0,\pi ],-2\le z\le 2.$$

The following recurrence relation for ${\mathsf{\Psi }}_{j\left(z\right)}$ is easily demonstrated:

$${\mathsf{\Psi }}_j\left(z\right)=z{\mathsf{\Psi }}_{j-1}\left(z\right)-{\mathsf{\Psi }}_{j-2}\left(z\right),j=2,3,\dots ,{\mathsf{\Psi }}_0\left(z\right)=2,{\mathsf{\Psi }}_1\left(z\right)=z.$$

A set of orthogonal polynomials defined on the interval $[0,1]$ is established by incorporating VLPs and $z=4\alpha -2$. These polynomials will be represented by the symbol ${\mathsf{\Phi }}_m\left(\beta \right)$, as illustrated in: ${\mathsf{\Phi }}_j\left(\alpha \right)=\mathsf{\Psi }(4\alpha -2)=\mathsf{\Psi }\left(2\sqrt{\alpha }\right).$

${\mathsf{\Phi }}_j\left(\alpha \right)$ have the following recurrence relation:

$${\mathsf{\Phi }}_{j+1}\left(\alpha \right)=(4\alpha -2){\mathsf{\Phi }}_{j-1}\left(\alpha \right)-{\mathsf{\Phi }}_{j-2}\left(\alpha \right),j=2,3,\dots ,$$

where, ${\mathsf{\Phi }}_0\left(\alpha \right)=2,{\mathsf{\Phi }}_1\left(\alpha \right)=4\alpha -2.$ Also, we find ${\mathsf{\Phi }}_j\left(0\right)=2{(-1)}^j$ and ${\mathsf{\Phi }}_j\left(1\right)=2,j=0,1,2,\dots .$

The analytical formula for ${\mathsf{\Phi }}_j\left(\alpha \right)$ is:

$${\mathsf{\Phi }}_k\left(\alpha \right)=2j\sum _{j=0}^k{(-1)}^j\frac{4^{k-j}\mathsf{\Gamma }(2k-j)}{\mathsf{\Gamma }(j+1)\mathsf{\Gamma }(2k-2j+1)}{\alpha }^{k-j},k=2,3,\dots .$$

The polynomials ${\mathsf{\Phi }}_k\left(\alpha \right)$ are orthogonal polynomials on $[0,1]$ w.r.t. $\frac{1}{\sqrt{\alpha -{\alpha }^2}}$, and so we have:

$$〈{\mathsf{\Phi }}_k\left(\alpha \right),{\mathsf{\Phi }}_j\left(\alpha \right)〉={\int }_0^1\frac{{\mathsf{\Phi }}_k\left(\alpha \right){\mathsf{\Phi }}_j\left(\alpha \right)}{\sqrt{\alpha -{\alpha }^2}}d\alpha =\left\{\begin{array}{cc}0,\hfill & k\ne j\ne 0,\hfill \\ 4\pi ,\hfill & k=j=0,\hfill \\ 2\pi ,\hfill & k=j\ne 0.\hfill \end{array}\right.$$

Let $\mathsf{\Omega }\left(\alpha \right)\in L^2[0,1]$, then using ${\mathsf{\Phi }}_k\left(\alpha \right)$, we have:

$$\mathsf{\Omega }\left(\alpha \right)=\sum _{k=0}^{\infty }{c}_k{\mathsf{\Phi }}_k\left(\alpha \right).$$

(1)

To convert $\mathsf{\Omega }\left(\alpha \right)$ into terms of ${\mathsf{\Phi }}_m\left(\alpha \right)$, it is necessary to evaluate $c_k.$ This can be expressed by considering only the first $m+1$ terms of (1), we obtain that

$${\mathsf{\Omega }}_m\left(\alpha \right)=\sum _{k=0}^m{c}_k{\mathsf{\Phi }}_k\left(\alpha \right),$$

(2)

The values of $c_k$ for $k=0,1,2,\dots ,m$ can be computed from:

$$c_0=\frac{1}{4\pi }{\int }_0^1\frac{\varphi \left(\alpha \right){\mathsf{\Phi }}_0\left(\alpha \right)}{\sqrt{\alpha -{\alpha }^2}}d\alpha ,c_k=\frac{1}{2\pi }{\int }_0^1\frac{\varphi \left(\alpha \right){\mathsf{\Phi }}_k\left(\alpha \right)}{\sqrt{\alpha -{\alpha }^2}}d\alpha ,k=1,2,3,\cdots .$$

(3)

2.2. Convergence Analysis

Lemma 1.

Assuming that the function ${\mathsf{\Omega }}_m\left(\alpha \right)$ belongs to $L^2[0,1]$ with the weight function $\frac{1}{\sqrt{\alpha -{\alpha }^2}}$ and $|{\mathsf{\Omega }}^{″}\left(\alpha \right)|\le \epsilon$, for some constant ε, the series given by (2) uniformly converges to $\mathsf{\Omega }\left(\alpha \right)$ as $m\to \infty$. Furthermore, the estimates listed below are accomplished:

1. 

Equation (2) has a bounded series of coefficients, meaning

$$\left|c_k\right|\le \frac{\epsilon }{4k\left(k^2-1\right)},k>2.$$

2. 

The norm for the error estimate $(L^2[0,1]$-norm) can be articulated as follows:

$$∥\mathsf{\Omega }\left(\alpha \right)-{\mathsf{\Omega }}_m\left(\alpha \right)∥<\frac{L}{12\sqrt{m^3}}.$$

3. 

If ${\mathsf{\Omega }}^{\left(m\right)}\left(\alpha \right)\in C[0,1]$, then the absolute error bound is satisfied

$$∥\mathsf{\Omega }\left(\alpha \right)-{\mathsf{\Omega }}_m\left(\alpha \right)∥\le \frac{\mathsf{\Delta }{\mathsf{\Pi }}^{m+1}}{(m+1)!}\sqrt{\pi },\mathsf{\Delta }=\underset{\alpha \in [0,1]}{max}{\mathsf{\Omega }}^{(m+1)}\left(\alpha \right),\mathrm{and}\mathsf{\Pi }=max\left\{1-{\alpha }_0,{\alpha }_0\right\}.$$

For a more in-depth exploration of the convergence analysis and the use of polynomials to approximate (2), readers may refer [33].

2.3. Nonstandard Finite Difference Scheme Notations

We provide the following elucidation of the discrete first derivative:

$$\begin{array}{c}\hfill \frac{d\mathsf{\Omega }}{d\beta }\to \frac{{\mathsf{\Omega }}_{s+1}-{\psi }_1\left(h\right){\mathsf{\Omega }}_s}{{\psi }_2\left(h\right)},\end{array}$$

(4)

where ${\psi }_1\left(h\right)$ and ${\psi }_2\left(h\right)$ are functions in the step-size discretization $h=\mathsf{\Delta }\left(\beta \right)$ and

$${\psi }_1\left(h\right)=1+O\left(h^2\right)\mathrm{and}{\psi }_2\left(h\right)=h+O\left(h^2\right).$$

The term used to denote the first derivative in the nonstandard finite difference technique is referred to as the presentation formula. Moreover, it is required that the denominator function adheres to the condition that $0<{\psi }_1\left(h\right)<1,h\to 0$. While the exact base function ${\psi }_2\left(h\right)$ is not specified, commonly used functions in the non-standard finite difference method can be introduced as follows: ${\psi }_2\left(h\right)=exp\left(h\right)-1,{\psi }_2\left(h\right)=h,{\psi }_2\left(h\right)=sinh\left(h\right),{\psi }_2\left(h\right)=\frac{1-exp(-\lambda h)}{\lambda },\mathrm{etc}.$ (see [34]).

3. Spectral Method for Solving CMSKDV

This section is subdivided into three sections, each delineating the fractional derivatives associated with distinct kernels. Recent works have investigated a fractional derivative based on the exponential decay rule, which can be described as a generalized power law function. A more generalized version of the exponential function called the Mittag-Leffler function served as the foundation for the fractional derivative with a non-local kernel that Dumitru Baleanu and Abdon Atangana presented. The depiction of intricate physical issues that follow the rules of power and exponential decay is made possible by this derivative [35,36].

The coupled KdV equations in the classical case are given as follows:

$$\begin{array}{c}\hfill {\varphi }_{\beta }-c{\varphi }_{\alpha \alpha \alpha }-{\mu }_2\varphi {\varphi }_{\alpha }-{\mu }_1\phi {\phi }_{\alpha }=0,\end{array}$$

(5)

and

$$\begin{array}{c}\hfill {\phi }_{\beta }-c{\phi }_{\alpha \alpha \alpha }+{\mu }_1\varphi {\phi }_{\alpha }=0,\end{array}$$

(6)

where $c,{\mu }_1$ and ${\mu }_2$ are constant parameters.

3.1. CMSKDV via the Liouville-Caputo Fractional Derivative

When replacing the classical derivative of the coupled KdV equations with a fractional derivative employing a power-law kernel or Liouville-Caputo fractional derivative (LC), the resulting system of equations for the coupled KdV reads as follows:

$$\begin{array}{c}\hfill {\varphi }_{\beta }-c{}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_1}\varphi -{\mu }_2\varphi {}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_2}\varphi -{\mu }_1\phi {}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_2}\phi =0\end{array}$$

(7)

and

$$\begin{array}{c}\hfill {\phi }_{\beta }-c{}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_1}\phi +{\mu }_1\varphi {}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_2}\phi =0,\end{array}$$

(8)

$$2<{\nu }_1\le 3,0<{\nu }_2\le 1,$$

where ${}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_1}$ and ${}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_2}$ are the LC fractional derivatives of order ${\nu }_1$ and ${\nu }_2$, respectively for the functions $\varphi \left(\alpha \right)$ and $\phi \left(\alpha \right)$ that belong to $H^1(0,b)$ defined in the following form [37,38].

$$\begin{array}{c}\hfill {}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_1}\phi \left(\alpha \right)=\frac{1}{\mathsf{\Gamma }(3-{\nu }_1)}{\int }_0^{\alpha }\frac{{\phi }^{\left(3\right)}\left(\tau \right)}{{(\alpha -\tau )}^{{\nu }_1-2}}d\tau (\alpha >0),\end{array}$$

(9)

$$\begin{array}{c}\hfill {}^{\mathrm{LC}}{D}_{\alpha }^{{\nu }_2}\phi \left(\alpha \right)=\frac{1}{\mathsf{\Gamma }(1-{\nu }_2)}{\int }_0^{\alpha }\frac{{\phi }^{\left(1\right)}\left(\tau \right)}{{(\alpha -\tau )}^{{\nu }_2}}d\tau (\alpha >0).\end{array}$$

(10)

The following theorem provides an approximate formula for the fractional derivative ${}^{\mathrm{LC}}{D}_{\alpha }^{\nu }\mathsf{\Omega }\left(\alpha \right),n-1<\nu \le n,$ with the power-law kernel.

Theorem 1

(see [39]). The LC-fractional-order derivative of the function ${\mathsf{\Omega }}_m\left(\alpha \right)$, as defined in Equation (2), can be determined using the following approximate formula:

$$\begin{array}{c}\hfill {}^{\mathrm{LC}}{D}_{\alpha }^{\nu }\left({\mathsf{\Omega }}_m\left(\alpha \right)\right)=\sum _{k=\lceil \nu \rceil }^m\sum _{𝓁=0}^{k-\mu }{\alpha }_k{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\nu ,𝓁}{}^{\mathrm{LC}}\left(\alpha \right),\end{array}$$

(11)

where

$$\begin{array}{c}\hfill {\vartheta }_{k,𝓁}=\frac{{(-1)}^{𝓁}{4}^{k-𝓁}2k\mathsf{\Gamma }(2k-𝓁)}{\mathsf{\Gamma }(𝓁+1)\mathsf{\Gamma }(2k-2𝓁+1)},\end{array}$$

(12)

and

$$\begin{array}{c}\hfill {\mathsf{\Pi }}_{\nu ,𝓁}{}^{\mathrm{LC}}\left(\alpha \right)=\frac{\mathsf{\Gamma }(k-𝓁+1)}{\mathsf{\Gamma }(k-𝓁-\nu +1)}{\alpha }^{k-𝓁-\nu }.\end{array}$$

(13)

Subsequently, we move forward with solving Equations (7) to (8) using the Vieta-Lucas spectral collocation technique, which is explained as follows:

$$\begin{array}{c}\hfill {\varphi }_m(\alpha ,\beta )=\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\end{array}$$

(14)

and

$$\begin{array}{c}\hfill {\phi }_m(\alpha ,\beta )=\sum _{k=0}^m{\phi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right).\end{array}$$

(15)

The two Equations (7) through (8) can be represented as follows in accordance with Theorem 1 and the expansions (14)–(15)

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\varphi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_k\left(\alpha \right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\nu ,𝓁}{}^{\mathrm{LC}}\left(\alpha \right)+{\mu }_2\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left(\alpha \right)\right)+{\mu }_1\left(\sum _{k=0}^m{\phi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \end{array}$$

$$\begin{array}{c}\hfill ·\left(\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left(\alpha \right)\right)\end{array}$$

(16)

and

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\phi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_k\left(\alpha \right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\nu ,𝓁}{}^{\mathrm{LC}}\left(\alpha \right)-{\mu }_1\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil \nu \rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left(\alpha \right)\right).\hfill \end{array}$$

(17)

For the Equations (16) and (17), we will evaluate them at m points ${\alpha }_p(p=0,1,\cdots ,m-1)$ as outlined below:

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\varphi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_k\left({\alpha }_p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)+{\mu }_2\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)\right)+{\mu }_1\left(\sum _{k=0}^m{\phi }_k\left(\beta \right){\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \end{array}$$

$$\begin{array}{c}\hfill ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)\right),\end{array}$$

(18)

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\phi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_k\left({\alpha }_p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\nu ,𝓁}{}^{\mathrm{LC}}{\left({\alpha }_p\right)}_p-{\mu }_1\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)\right).\hfill \end{array}$$

(19)

Moreover, through substitution from Equation (2), we can articulate the corresponding boundary conditions, leading to the ensuing set of equations:

$$\begin{array}{c}\hfill \sum _{k=0}^m{(-1)}^k{\varphi }_k\left(\beta \right)=g_1\left(\beta \right),\sum _{k=0}^m{(-1)}^k{\phi }_k\left(\beta \right)=h_1\left(\beta \right),\end{array}$$

(20)

$$\begin{array}{c}\hfill \sum _{k=0}^m{\varphi }_k\left(\beta \right)=g_2\left(\beta \right),\sum _{k=0}^m{\phi }_k\left(\beta \right)=h_2\left(\beta \right)\end{array}$$

(21)

and

$$\begin{array}{c}\hfill \sum _{k=0}^{m}2k^2{\varphi }_k\left(\beta \right)=g_3\left(\beta \right)\mathrm{and}\sum _{k=0}^{m}2k^2{\phi }_k\left(\beta \right)=h_3\left(\beta \right).\end{array}$$

(22)

Ultimately, applying the nonstandard finite difference approach as described in [39] leads us to a set of nonlinear algebraic equations, we obtain

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{{\varphi }_k^p-{\varphi }_k^{p-1}}{0.5(exp(2h)-1)}{\mathsf{\Phi }}_k\left({\alpha }_p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\varphi }_k^p{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)+{\mu }_2\left(\sum _{k=0}^m{\varphi }_k^p{\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\varphi }_k^p{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)\right)+{\mu }_1\left(\sum _{k=0}^m{\phi }_k^p{\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k^p{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)\right)=0,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{{\phi }_k^p-{\phi }_k^{p-1}}{0.5(exp(2h)-1)}{\mathsf{\Phi }}_k\left({\alpha }_p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k^p{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\nu ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)-{\mu }_1\left(\sum _{k=0}^m{\varphi }_k^p{\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k^p{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{\theta ,𝓁}{}^{\mathrm{LC}}\left({\alpha }_p\right)\right).\hfill \end{array}$$

(24)

These can be solved using well-known methods such as the Newton-Raphson technique, enabling us to ascertain the values of ${\varphi }_k^p$ and ${\phi }_k^p$ for $k=0,1,\cdots ,m.$

3.2. CMSKDV via the Caputo-Fabrizio Fractional Derivative

By replacing the classical derivative with the fractional derivative incorporating either the exponential-decay kernel or the Caputo-Fabrizio (CF),the coupled (KdV) equations can be represented in the following manner:

$$\begin{array}{c}\hfill {\varphi }_{\beta }-c{}^{\mathrm{CF}}{D}_{\alpha }^{{\nu }_1}\varphi -{\mu }_2\varphi {}^{\mathrm{CF}}{D}_{\alpha }^{{\nu }_2}\varphi -{\mu }_1\phi {}^{\mathrm{CF}}{D}_{\alpha }^{{\nu }_2}\phi =0\end{array}$$

(25)

and

$$\begin{array}{c}\hfill {\phi }_{\beta }-c{}^{\mathrm{CF}}{D}_{\alpha }^{{\nu }_1}\phi +{\mu }_1\varphi {}^{\mathrm{CF}}{D}_{\alpha }^{{\nu }_2}\phi =0.\end{array}$$

(26)

The fractional derivative operator ${}^{\mathrm{CF}}{D}_{\alpha }^{\nu }$ ($\nu$ might be ${\nu }_1$ or ${\nu }_2$) of order $n<\nu <n+1$, which incorporates the exponential-decay kernel, is defined as follows [40]:

$$\begin{array}{c}\hfill {}^{\mathrm{CF}}{D}_{a+}^{\nu }\phi \left(\alpha \right)={}^{\mathrm{CF}}{D}_{a+}^{\kappa }\left(D^n\phi \left(\alpha \right)\right)=\frac{M\left(\kappa \right)}{1-\kappa }{\int }_a^{\alpha }{\phi }^{(n+1)}\left(\tau \right)e^{-\frac{\kappa (\alpha -\tau )}{1-\kappa }}d\tau ,\end{array}$$

(27)

where $n=\lfloor \nu \rfloor$ = the floor of $\nu$ (that is, the integer part) and $\kappa =\lceil \nu \rceil$ = the ceiling of $\nu$ (that is, the decimal part) and $M\left(\kappa \right)$ is the normalization function.

Theorem 2

(see [30,40]). The left-hand side features a fractional derivative employing an exponential-decay kernel with an order $\nu \in (n,n+1)$

$$\phi \left(\alpha \right)={\alpha }^p(p\geqq \lceil \nu \rceil )$$

can be represented in the following manner:

$$\begin{array}{cc}\hfill & {}^{\mathrm{CF}}{D}_{0+}^{\nu }{\alpha }^p=\frac{M\left(\kappa \right)\mathsf{\Gamma }(p+1)}{1-\kappa }\left(\sum _{i=0}^{p-n-1}\frac{{(-1)}^i{\alpha }^{p-n-1-i}}{\mathsf{\Gamma }(p-n-i){\left(\frac{\kappa }{1-\kappa }\right)}^{i+1}}+\frac{{(-1)}^{p-n}}{{\left(\frac{\kappa }{1-\kappa }\right)}^{p-n}}{e}^{-\frac{\kappa \alpha }{1-\kappa }}\right).\hfill \end{array}$$

(28)

The formula in Theorem 1 of [40] when $a=0$ corresponds to a specific case of the Formula (28). For additional details regarding the definitions and characteristics of fractional CF-derivatives, refer to [40].

The following theorem provides the approximation formula of ${}^{\mathrm{CF}}{D}_{0+}^{\nu }{\mathsf{\Omega }}_m\left(x\right)$ by means of (1).

Theorem 3

(see [30]). The ${}^{\mathrm{CF}}{D}_{\alpha }^{\nu }\left({\mathsf{\Omega }}_m\left(\alpha \right)\right)$ can be written as follows:

$$\begin{array}{c}\hfill {}^{\mathrm{CF}}{D}_{0+}^{\nu }\left({\mathsf{\Omega }}_m\left(\alpha \right)\right)=\sum _{k=\lceil \nu \rceil }^m\sum _{𝓁=\lceil \nu \rceil }^k{a}_k{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left(\alpha \right),\end{array}$$

(29)

where

$$\begin{array}{c}\hfill {\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left(\alpha \right)=\left(\frac{M\left(\kappa \right)}{1-\kappa }\right)\left(\frac{{(-1)}^{𝓁-n}}{(\frac{\kappa }{1-\kappa }{)}^{𝓁-n}}{e}^{-\frac{\kappa \alpha }{1-\kappa }}+\sum _{p=0}^{𝓁-n-1}\frac{{(-1)}^p{\alpha }^{𝓁-n-1-p}}{\mathsf{\Gamma }(𝓁-n-p)(\frac{\kappa }{1-\kappa }{)}^{p+1}}\right).\end{array}$$

(30)

By using Equations (25) and (26), (29) and the Formula (2) with CF, we obtain

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\varphi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_k\left(\alpha \right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left(\alpha \right)+{\mu }_2\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left(\alpha \right)\right)+{\mu }_1\left(\sum _{k=0}^m{\phi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left(\alpha \right){\alpha }^{k-𝓁-\theta }\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\phi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_j\left(\alpha \right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k\left(\beta \right),{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left(\alpha \right)-{\mu }_1\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left(\alpha \right)\right).\hfill \end{array}$$

(32)

To obtain the ultimate expressions denoted as (31) and (32), we will assess them at m specific points ${\alpha }_p;(p=0,1,\cdots ,m-1)$ according to the following procedure:

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{{\varphi }_k^p-{\varphi }_k^{p-1}}{0.5(exp(2h)-1)}{\mathsf{\Phi }}_k\left({\alpha }_p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\varphi }_k\left(\beta \right),{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left({\alpha }_p\right)+{\mu }_2\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left({\alpha }_p\right)\right)+{\mu }_1\left(\sum _{k=0}^m{\phi }_k\left(\beta \right){\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left({\alpha }_p\right)\right),\hfill \end{array}$$

(33)

and

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{{\phi }_k^p-{\phi }_k^{p-1}}{0.5(exp(2h)-1)}{\mathsf{\Phi }}_j\left({\alpha }_p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left({\alpha }_p\right)-{\mu }_1\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ \hfill & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{CF}}\left({\alpha }_p\right)\right).\hfill \end{array}$$

(34)

In Equations (33) to (34),the nonstandard finite difference approach as described in [39] was utilized to compute the derivative concerning with $\beta$. Consequently, the equations were transformed into a set of algebraic equations.

One can solve (33) to (34) using established methods like the Newton-Raphson technique, allowing us to determine the values of ${\varphi }_k^p$ and ${\phi }_k^p$ for $k=0,1,\cdots ,m.$

3.3. CMSKDV via the Atangana-Baleanu-Caputo Fractional Derivative

When we replace the classical derivative with the fractional derivative incorporating either the Mittag-Leffler kernel or the Atangana-Baleanu-Caputo fractional-order derivative (ABC), the resultant coupled Korteweg-de Vries system of equations is represented as follows:

$$\begin{array}{c}\hfill {\varphi }_{\beta }-c{}^{\mathrm{ABC}}{D}_{\alpha }^{{\nu }_1}\varphi -{\mu }_2\varphi {}^{\mathrm{ABC}}{D}_{\alpha }^{{\nu }_2}\varphi -{\mu }_1\phi {}^{\mathrm{ABC}}{D}_{\alpha }^{{\nu }_2}\phi =0\end{array}$$

(35)

and

$$\begin{array}{c}\hfill {\phi }_{\beta }-c{}^{\mathrm{ABC}}{D}_{\alpha }^{{\nu }_1}\phi +{\mu }_1\varphi {}^{\mathrm{ABC}}{D}_{\alpha }^{{\nu }_2}\phi =0.\end{array}$$

(36)

The fractional derivative ${{}_a^{\mathrm{ABC}}D}^{\nu }$ ($\nu$ can be ${\nu }_1$ or ${\nu }_2$.) of order $n<\nu \leqq n+1$ is defined by [36,41].

$$\begin{array}{cc}\hfill {{}_0^{\mathrm{ABC}}D}^{\nu }\phi \left(\alpha \right)& ={{}_0^{\mathrm{ABC}}D}^{\kappa }\left(D^n\phi \left(\alpha \right)\right)\hfill \\ \hfill & =\frac{M\left(\kappa \right)}{1-\kappa }{\int }_0^{\alpha }{\phi }^{(n+1)}\left(\tau \right)E_{\kappa }\left(-\frac{\kappa }{1-\kappa }{(\alpha -\tau )}^{\kappa }\right)d\tau ,\hfill \end{array}$$

where $n=\lfloor \nu \rfloor$ = the floor of $\nu$ (that is, the integer part) and $\kappa =\lceil \nu \rceil$ = the ceiling of $\nu$ (that is, the decimal part).

$M\left(\kappa \right)$ is the normalization function such that $M\left(0\right)=M\left(1\right)=1$ and

$$E_{\nu }\left(\alpha \right)=\sum _{k=0}^{\infty }\frac{{\alpha }^k}{\mathsf{\Gamma }(k\nu +1)}$$

is the Mittag-Leffler function [42] Section 4.

Theorem 4

(see [31]). The ABC-derivative of order $\nu \in (n,n+1]$ of the following function:

$$\phi \left(\alpha \right)={\alpha }^p(p\geqq \lceil \nu \rceil )$$

is given by

$$\begin{array}{c}\hfill {{}_0^{\mathrm{ABC}}D}^{\nu }{\alpha }^p=\frac{M\left(\kappa \right)\mathsf{\Gamma }(p+1)}{1-\kappa }\sum _{i=0}^{\infty }{\left(-\frac{\kappa }{1-\kappa }\right)}^i\frac{{\alpha }^{p+i\kappa -n}}{\mathsf{\Gamma }(p+i\kappa -n+1)}.\end{array}$$

(37)

In the following theorem, we present the fundamental approximation formula for ${{}_0^{\mathrm{ABC}}D}^{\nu }{\mathsf{\Omega }}_m\left(\alpha \right)$ by utilizing the approximation (2).

Theorem 5

(see [31]). The ABC-derivative ${{}_0^{\mathrm{ABC}}D}^{\nu }\left({\mathsf{\Omega }}_m\left(\alpha \right)\right)$ can be expressed as follows:

$${{}_0^{\mathrm{ABC}}D}^{\nu }\left({\mathsf{\Omega }}_m\left(\alpha \right)\right)=\sum _{k=\lceil \nu \rceil }^m\sum _{𝓁=\lceil \nu \rceil }^{k-\lceil \nu \rceil }{\alpha }_k{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left(\alpha \right),$$

(38)

where

$$\begin{array}{cc}\hfill & {\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left(\alpha \right)=\frac{M\left(\kappa \right)\mathsf{\Gamma }(𝓁+1)}{1-\kappa }\sum _{s=0}^{\infty }\frac{{\left(-\frac{\kappa }{1-\kappa }\right)}^s{\alpha }^{𝓁+p\kappa -n}}{\mathsf{\Gamma }(𝓁+s\kappa -n+1)}.\hfill \end{array}$$

(39)

It can be readily noticed that

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\varphi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_k\left(\alpha \right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left(\alpha \right)+{\mu }_2\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left(\alpha \right)\right)+{\mu }_1\left(\sum _{k=0}^m{\phi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left(\alpha \right){\alpha }^{k-𝓁-\theta }\right),\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{d{\phi }_k\left(\beta \right)}{d\beta }{\mathsf{\Phi }}_k\left(\alpha \right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k\left(\beta \right),{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left(\alpha \right)-{\mu }_1\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \end{array}$$

$$\begin{array}{c}\hfill ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left(\alpha \right)\right).\end{array}$$

(41)

To derive the final formulations indicated as (40) and (41), we will evaluate them at m distinct points ${\alpha }_p;(p=0,1,\cdots ,m-1)$ following the outlined procedure:

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{{\varphi }_k^p-{\varphi }_k^{p-1}}{0.5(exp(2h)-1)}{\mathsf{\Phi }}_k\left({\alpha }_p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\varphi }_k\left(\beta \right),{\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left({\alpha }_p\right)+{\mu }_2\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left(\alpha \right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\varphi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left({\alpha }_p\right)\right)+{\mu }_1\left(\sum _{k=0}^m{\phi }_k\left(\beta \right){\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left({\alpha }_p\right)\right),\hfill \end{array}$$

(42)

and

$$\begin{array}{cc}\hfill \sum _{k=0}^m\frac{{\varphi }_k^p-{\varphi }_k^{p-1}}{0.5(exp(2h)-1)}{\mathsf{\Phi }}_k\left(\alpha p\right)& =c\sum _{k=\lceil {\nu }_1\rceil }^m\sum _{𝓁=\lceil {\nu }_1\rceil }^{k-\lceil {\nu }_1\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left({\alpha }_p\right)-{\mu }_1\left(\sum _{k=0}^m{\varphi }_k\left(\beta \right){\mathsf{\Phi }}_k\left({\alpha }_p\right)\right)\hfill \\ & ·\left(\sum _{k=\lceil {\nu }_2\rceil }^m\sum _{𝓁=\lceil {\nu }_2\rceil }^{k-\lceil {\nu }_2\rceil }{\phi }_k\left(\beta \right){\vartheta }_{k,𝓁}{\mathsf{\Pi }}_{k,𝓁,\kappa }{}^{\mathrm{ABC}}\left({\alpha }_p\right)\right).\hfill \end{array}$$

(43)

These equations can be resolved utilizing well-known approaches such as the Newton-Raphson technique, enabling us to calculate the values of ${\varphi }_k^p$ and ${\phi }_k^p$ for $k=0,1,\cdots ,m.$

4. Numerical Results and Discussion

In this section, we utilize the three previously mentioned approaches to perform a numerical evaluation of the multi-space fractional-order coupled Korteweg-de Vries equation, incorporating three different kernels. Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 are included to visually represent the numerical results. The numerical results provided in this section will be computed using the designated values: $c=-1,$ ${\mu }_1=-6,$ ${\mu }_2=7,$ and $a=1.$

In Figure 1, a comparison is made between the exact solution and the approximate solution of (7) and (8) in the case of the LC-derivative, as specified above with ${\nu }_1=2.8,$ and ${\nu }_2=0.8.$

The exact solutions of the coupled Korteweg-de Vries equation under classical conditions are given as follows:

$$\begin{array}{c}\hfill \varphi (\alpha ,\beta )=\phi (\alpha ,\beta )=\frac{4c^{2}exp\left(c\left(\alpha -\beta c^2\right)\right)}{{\left(exp\left(c\left(\alpha -\beta c^2\right)\right)+1\right)}^2}.\end{array}$$

(44)

In order to apply the initial condition in the exact solution for the coupled Korteweg-de Vries equation, we set $\beta =0$. Subsequently, we acquire the initial values as follows:

$$\begin{array}{c}\hfill \varphi \left(\alpha \right)=\phi \left(\alpha \right)=\frac{4c^{2}exp\left(c\left(\alpha \right)\right)}{{\left(exp\left(c\left(\alpha \right)\right)+1\right)}^2}\end{array}$$

(45)

Additionally, the boundary conditions, once $\alpha =0$ and $\alpha =1$ are set, are expressed as follows:

$$\begin{array}{c}\hfill \varphi (0,\beta )=\phi (0,\beta )=\frac{4c^{2}exp\left(c\left(-\beta c^2\right)\right)}{{\left(exp\left(c\left(-\beta c^2\right)\right)+1\right)}^2},\end{array}$$

(46)

$$\begin{array}{c}\hfill \varphi (1,\beta )=\phi (1,\beta )=\frac{4c^{2}exp\left(c\left(1-\beta c^2\right)\right)}{{\left(exp\left(c\left(1-\beta c^2\right)\right)+1\right)}^2},\end{array}$$

(47)

$$\begin{array}{c}\hfill {\varphi }_{\alpha }(0,\beta )={\phi }_{\alpha }(0,\beta )=\frac{4c^3{e}^{c\left(\beta \left(-c^2\right)\right)}}{{\left(e^{c\left(\beta \left(-c^2\right)\right)}+1\right)}^2}-\frac{8c^3{e}^{2c\left(\beta \left(-c^2\right)\right)}}{{\left(e^{c\left(\beta \left(-c^2\right)\right)}+1\right)}^3.}\end{array}$$

(48)

In all calculations conducted in this study, the initial and boundary conditions mentioned earlier will be utilized.

In Figure 2, we display the absolute discrepancy between the precise solution (44) and the numerical solution (7) and (8) for the identical set of parameters as depicted in Figure 1.

Given the absence of an exact solution for the fractional-order equation in this scenario, we will assess the absolute error between two consecutive steps.

Figure 3 illustrates the plot of the discrepancy between two consecutive stages.

Based on the data presented in the preceding three figures, it is evident that the numerical solutions we have presented exhibit high accuracy with a remarkably low error rate. Our emphasis on fractional-order equations stems from their effectiveness in gauging the error between two consecutive phases, aligning with the primary objective of our study.

For the LC-derivative, we have carried out the necessary calculations and generated the corresponding data representations.

We conducted a comparable investigation to that of LC, as depicted in Figure 4, Figure 5 and Figure 6, but this time incorporating the CF-derivative operator.

We have also examined the numerical solutions of Equations (42)–(43) in Figure 7, Figure 8 and Figure 9, this time employing the ABC-derivative operator.

Based on this analysis, we have noted a relatively minimal numerical calculation error. As a result, this study’s findings can be extended and applied to other fractional-order equations and systems.

In

Table 1. The absolute error between Equations (7) and (8) and the exact solution (44) is represented by the first column and the second column while the absolute error of two steps of (7) and (8) is represented by the third column and the fourth column for ${\nu }_1=2.8,{\nu }_2=0.8,c=1,{\mu }_1=7,{\mu }_2=-6,$ and $a=1.$

$\mathit{\alpha }$	$|\mathit{\varphi }(\mathit{\alpha },\mathit{\beta })-{\mathit{\varphi }}_{\mathit{LC}}({\mathit{\alpha }}_{\mathit{n}},{\mathit{\beta }}_{\mathit{n}})|$	$|\mathit{\phi }(\mathit{\alpha },\mathit{\beta })-{\mathit{\phi }}_{\mathit{LC}}({\mathit{\alpha }}_{\mathit{n}},{\mathit{\beta }}_{\mathit{n}})|$	$|{\mathit{\varphi }}_{\mathit{LC},\mathit{n}}-{\mathit{\varphi }}_{\mathit{LC},\mathit{n}+1}|$	$|{\mathit{\phi }}_{\mathit{LC},\mathit{n}}-{\mathit{\phi }}_{\mathit{LC},\mathit{n}+1}|$
0	$1.11022\times {10}^{-16}$	$1.11022\times {10}^{-16}$	$4.97499\times {10}^{-11}$	$4.97495\times {10}^{-11}$
0.2	$8.30951\times {10}^{-5}$	$2.35526\times {10}^{-5}$	$1.12725\times {10}^{-6}$	$4.48712\times {10}^{-7}$
0.4	$1.35067\times {10}^{-4}$	$4.38737\times {10}^{-5}$	$1.69840\times {10}^{-6}$	$7.35767\times {10}^{-7}$
0.6	$2.28422\times {10}^{-4}$	$9.50490\times {10}^{-5}$	$2.41060\times {10}^{-6}$	$1.11918\times {10}^{-6}$
0.8	$3.85878\times {10}^{-4}$	$1.56810\times {10}^{-4}$	$4.12192\times {10}^{-6}$	$1.81449\times {10}^{-6}$
1	$1.11022\times {10}^{-16}$	$1.11022\times {10}^{-16}$	$3.63416\times {10}^{-7}$	$3.63416\times {10}^{-7}$

,

Table 2. In the case of the CF-derivative, the absolute error between Equations (25) and (26) and the exact solution (44) is represented by the first column and the second column while the absolute error of two steps of (25) and (26) is represented by the third column and the fourth column for ${\nu }_1=2.8,{\nu }_2=0.8,c=1,{\mu }_1=7,{\mu }_2=-6,$ and $a=1.$

$\mathit{\alpha }$	$|\mathit{\varphi }(\mathit{\alpha },\mathit{\beta })-{\mathit{\varphi }}_{\mathit{CF}}({\mathit{\alpha }}_{\mathit{n}},{\mathit{\beta }}_{\mathit{n}})|$	$|\mathit{\phi }(\mathit{\alpha },\mathit{\beta })-{\mathit{\phi }}_{\mathit{CF}}({\mathit{\alpha }}_{\mathit{n}},{\mathit{\beta }}_{\mathit{n}})|$	$|{\mathit{\varphi }}_{\mathit{CF},\mathit{n}}-{\mathit{\varphi }}_{\mathit{CF},\mathit{n}+1}|$	$|{\mathit{\phi }}_{\mathit{CF},\mathit{n}}-{\mathit{\phi }}_{\mathit{CF},\mathit{n}+1}|$
0	$1.92593\times {10}^{-34}$	0	$4.97498\times {10}^{-11}$	$4.97501\times {10}^{-11}$
0.2	$2.75570\times {10}^{-5}$	$1.59500\times {10}^{-4}$	$3.27900\times {10}^{-7}$	$1.07397\times {10}^{-6}$
0.4	$1.89708\times {10}^{-4}$	$4.38737\times {10}^{-5}$	$2.55455\times {10}^{-7}$	$1.27453\times {10}^{-6}$
0.6	$2.49232\times {10}^{-6}$	$1.04730\times {10}^{-4}$	$1.94268\times {10}^{-7}$	$8.70804\times {10}^{-7}$
0.8	$8.67747\times {10}^{-6}$	$1.73810\times {10}^{-4}$	$2.06507\times {10}^{-7}$	$1.54168\times {10}^{-6}$
1	$3.851860\times {10}^{-34}$	$1.11022\times {10}^{-16}$	$3.63417\times {10}^{-7}$	$3.63417\times {10}^{-7}$

and

Table 3. In the case of the ABC-derivative, the absolute error between Equations (40) and (41) and the exact solution (44) is represented by the first column and the second column while the absolute error of two steps of (40) and (41) is represented by the third column and the fourth column for ${\nu }_1=2.8,{\nu }_2=0.8,c=1,{\mu }_1=7,{\mu }_2=-6,$ and $a=1.$

$\mathit{\alpha }$	$|\mathit{\varphi }(\mathit{\alpha },\mathit{\beta })-{\mathit{\varphi }}_{\mathit{ABC}}({\mathit{\alpha }}_{\mathit{n}},{\mathit{\beta }}_{\mathit{n}})|$	$|\mathit{\phi }(\mathit{\alpha },\mathit{\beta })-{\mathit{\phi }}_{\mathit{ABC}}({\mathit{\alpha }}_{\mathit{n}},{\mathit{\beta }}_{\mathit{n}})|$	$|{\mathit{\varphi }}_{\mathit{CF},\mathit{n}}-{\mathit{\varphi }}_{\mathit{ABC},\mathit{n}+1}|$	$|{\mathit{\phi }}_{\mathit{CF},\mathit{n}}-{\mathit{\phi }}_{\mathit{ABC},\mathit{n}+1}|$
0	$1.11022\times {10}^{-16}$	0	$4.97496\times {10}^{-11}$	$4.97498\times {10}^{-11}$
0.2	$1.75119\times {10}^{-5}$	$1.50089\times {10}^{-4}$	$1.29273\times {10}^{-7}$	$1.11919\times {10}^{-6}$
0.4	$2.283092244170781\times {10}^{-5}$	$1.54883\times {10}^{-4}$	$5.88663\times {10}^{-8}$	$1.09897\times {10}^{-6}$
0.6	$1.41835\times {10}^{-6}$	$3.42845\times {10}^{-6}$	$1.25380\times {10}^{-7}$	$1.04458\times {10}^{-7}$
0.8	$3.190648624473713\times {10}^{-5}$	$1.29585\times {10}^{-5}$	$3.09597\times {10}^{-8}$	$1.75898\times {10}^{-7}$
1	0	0	$3.63415\times {10}^{-7}$	$3.63417\times {10}^{-7}$

, the tables illustrate the absolute error between the approximate solutions of two equations under the LC-derivative, CF-derivative, and ABC-derivative. Additionally, the tables provide insights into the extent of absolute error between consecutive steps when encountering a non-integer order for the LC-derivative, CF-derivative, and ABC-derivative. The results from these tables indicate that the order of the error is notably minimal. This observation serves as a strong indication of the precision and efficacy of the approach presented in this study.

5. Conclusions

The problem of multi-space fractional-order coupled Korteweg-de Vries, which involves various kernels, was successfully converted into a system of differential equations. Subsequently, these equations were converted into algebraic forms employing the non-standard finite difference approach.

We employed the commonly used numerical technique, often referred to as the Newton-Raphson method, to address the resultant set of nonlinear algebraic equations. The precision of the examined approximations was confirmed by computing the absolute error between the exact solutions and the approximations. Additionally, we determined the absolute error between successive approximations in the scenario of non-integer orders.

We verified the accuracy of the analyzed approximations by calculating the absolute discrepancy between the exact solutions and the approximations. Furthermore, we assessed the absolute error between consecutive approximations in cases involving non-integer orders.

In future work, the current study’s findings can be extended to analyze the behavior of positive solutions [43,44]. The scope of the investigation may also be broadened to encompass space-time. Furthermore, the study discussed in this paper can be extended to cover other special functions, such as Legendre and Chebyshev polynomials. In this context, recent research has been conducted on the subject of fractional derivatives and special functions, as in [45,46,47,48].