Global attractivity of a rational difference equation with higher order and its applications

1.   Introduction

In various fields of science, nonlinear evolution equations practically model many natural, biological and engineering processes. For example, PDEs are very popular and are used in physics to study traveling wave solutions. They have played a crucial role in illustrating the nature of nonlinear problems. PDEs are collected to control the diffusion of chemical reactions. In biology, they play a fundamental role in describing various phenomena, such as population growth. In addition, natural phenomena such as fluid dynamics, plasma physics, optics and optical fibers, electromagnetism, quantum mechanics, ocean waves, and others are studied using PDEs. The qualitative and quantitative characteristics of these phenomena can be identified from the behaviors and shapes of their solutions. Therefore, finding the analytic solutions to such phenomena is a fundamental topic in mathematics. Scientists have developed sparse fundamental approaches to find analytic solutions for nonlinear PDEs. Among these techniques, I present integration methods from ^[1] and ^[2], the modified F-expansion and Generalized Algebraic methods, respectively. Bekir and Unsal ^[3] proposed the first integral method to find the analytical solution of nonlinear equations. Kumar, Seadawy and Joardar ^[4] used the improved Kudryashov technique to extract fractional differential equations. Adomian ^[5] proposed the Adomian decomposition technique to find the solution of frontier problems of physics. ^[6] uses an exploratory method to find explicit solutions of non-linear PDEs. Many different methods of solving equations arising from natural phenomena and some of their analytic solutions, such as dark and light solitons, non-local rogue waves, an occasional wave and mixed soliton solutions, are exhibited and can be found in ^{[7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40]}.

The Novikov-Veselov (NV) system ^[41,42] is given by

where $\alpha, \, \beta, \, \gamma$ and $\lambda$ are constants. Barman ^[42] declared that Eq (1.1) is involved to represent tidal and tsunami waves, electro-magnetic waves in communication cables and magneto-sound and ion waves in plasma. In ^[42], the generalized Kudryashov method was utilized to have traveling wave solutions for Eq (1.1). According to Croke ^[43], the Novikov-Veselov system is generalized from the KdV equations which were examined by Novikov and Veselov. Croke ^[43] used several approaches, (the extended mapping, the Hirota and the extended tanh-function approaches) in the proposed system to achieve numerous soliton solutions, such as breathers, and constrained analytic solutions. Boiti, Leon, and Manna ^[44] applied the inverse dispersion technique to solve (1.1) for a particular type of initial value. Numerical solutions and a study of the stability of solutions for the proposed equation were presented by Kazeykina and Klein ^[45]. The Nizhnik-Novikov-Veselov system for two dimensions was also solved using the Kansa technique to find the numerical results ^[46]. To the best of my knowledge, the stability and error analysis of the numerical scheme presented here has not yet been discussed for system (1.1). Therefore, this has motivated me enormously to do so. The primary purpose is to obtain multiple analytic solutions to system (1.1) by using both the modified F-expansion and Generalized Algebraic methods. In connection with the numerical solution, the method of finite differences is utilized to achieve numerical results for the studied system. I graphically and analytically compare the traveling wave solutions and numerical results. Undoubtedly, the presented results strongly contribute to describing physical problems in practice.

The outline of this article is provided in this paragraph. Section 2 summarizes the employed methods. All the analytic solutions are extracted in Section 3. The shooting and BVP results for the proposed system are presented in Section 4. In addition, I examine the numerical solution of the system (1.1) in Section 5. Sections 6 and 7 study the stability and error analysis of the numerical scheme, respectively. Section 8 presents the results and discussion.

2.   Summary of proposed methods

Considering the development equation with physical fields $\Psi(x, y, t)$, $\Phi(x, y, t)$ and $\Gamma(x, y, t)$ in the variables $x,$ $y$ and $t$ is given in the following form:

Step 1. We extract the traveling-wave solutions of System (1.1) that are formed as follows:

where $w$ is the wave speed.

Step 2. The nonlinear evolution (2.1) is reduced to the following ODE:

where $Q_2$ is a polynomial in $\psi(\eta), \; \phi(\eta),$ $\Theta(\eta)$ and their total derivatives.

2.1.   The modified F-expansion method

According to the modified F-expansion method, the solutions of (2.3) are given by the form

and $F(\eta)$ is a solution of the following differential equation:

where $\mu_0$, $\mu_1$, $\mu_2$, are given in Table 1 ^[1], and $\rho_k$, $q_k$ are to be determined later.

2.2.   The generalized algebraic method

According to the generalized direct algebraic method, the solutions of (2.3) are given by

and $G(\eta)$ is a solution of the following differential equation:

where $\nu_k,$ and $r_k$ are to be determined, and $N$ is an integer number obtained by the highest degree of the nonlinear terms and the highest order of the derivatives. $\varepsilon$ is user-specified, usually taken with $\varepsilon = \pm1$, and $\delta_k$, $k = 0, 1, 2, 3, 4,$ are given in Table 2 ^[2].

3.   Methodology

Consider the Novikov-Veselov (NV) system

a system of PDEs in the unknown functions $\Psi = \Psi(x, y, t), \; \Phi = \Phi(x, y, t)$, $\Gamma = \Gamma(x, y, t)$ and their partial derivatives. I plug the transformations

into Eq (3.1) to reduce it to a system of ODEs given by

Integrating $\Theta_\eta = \psi_\eta$ and $\phi_\eta = \psi_\eta$ yields

Substituting (3.4) into the first equation of (3.3) and integrating once with respect to $\eta$ yields

Balancing $\psi_{\eta\eta}$ with $\psi^2$ in (3.5) calculates the value of $N = 2.$

3.1.   Application of modified F-expansion method

According to the modified F-expansion method with $N = 2$, the solutions of (3.5) are

and $F(\eta)$ is a solution of the following differential equation:

where $\mu_0$, $\mu_1$, $\mu_2$ are given in Table 1. To explore the analytic solutions to (3.5), I ought to follow the subsequent steps.

$\textbf{Step 1.}$ Placing (3.6) along with (3.7) into Eq (3.5) and gathering the coefficients of $F(\eta)^j$, $j$ = $-4, -3, -2, -1, 0, 1, 2, 3, 4,$ to zeros gives a system of equations for $\rho_0, \; \rho_k, \; q_k$, $k = 1, 2.$

$\textbf{Step 2.}$ Solve the resulting system using mathematical software: for example, Mathematica or Maple.

$\textbf{Step 3.}$ Choosing the values of $\mu_0, \; \mu_1$ and $\mu_2$ and the function $F(\eta)$ from Table 1 and substituting them along with $\rho_0, \; \rho_k, \; q_k$, $k = 1, 2,$ in (3.6) produces a set of trigonometric function and rational solutions to (3.5).

Applying the above steps, I determine the values of $\rho_0, \; \rho_1, \; \rho_2, \; q_1, \; q_2$ and $w$ as follows:

$(1).$ When $\mu_0 = 0,$ $\mu_1 = 1$ and $\mu_2 = -1$, I have two cases.

$\textbf{Case 1.}$

The solution is given by

$\textbf{Case 2.}$

The solution is given by

Figure 1 presents the time evolution of the analytic solutions $(a)$ $\Psi_1$ and $(b)$ $\Psi_2$ with $t = 0, 10, 20$. The parameter values are $x0 = -20, \ \alpha = 0.50, \beta = 0.6, \ \gamma = -1.5,$ and $\lambda = 1.$ Figure 2 presents the wave behavior by changing a certain parameter value and fixing the values of the others. Figure 2$(a, b)$ presents the behavior of $\Psi_1$ when I change the values of (a) $\alpha$ or $\beta$ and (b) $\gamma$ or $\lambda$. In Figure 2$(a)$ it can also be seen that the value of $\alpha$ or $\beta$ affects the direction and amplitude of the wave, such that a negative value always makes the wave negative, its amplitude decreases when $\alpha, \, \beta \to 0$, and its amplitude increases when $\alpha, \, \beta \to \infty.$ In Figure 2$(b)$ the value of $\gamma$ or $\lambda$ affects the direction and amplitude of the wave, such that a negative value always makes the wave negative, and its amplitude decreases when the value of $\gamma$ or $\lambda$ increases. In Figure 2, $(c)$ and $(d)$ present the wave behavior of $\Psi_2.$

$(2).$ When $\mu_0 = 0,$ $\mu_1 = -1,$ and $\mu_2 = 1$, I have two cases.

$\textbf{Case 3.}$ The solution is given by

$\textbf{Case 4.}$ The solution is given by

$(3).$ When $\mu_0 = 1,$ $\mu_1 = 0,$ and $\mu_2 = -1$, I have

$\textbf{Case 5.}$ The solution is given by

$(4).$ When $\mu_0 = \pm 1,$ $\mu_1 = 0,$ and $\mu_2 = \pm 1$, I have one case.

$\textbf{Case 6.}$ The solution is given by

3.2.   Application of the generalized algebraic method

According to the generalized algebraic method, the solutions of (3.5) are given by the form

where $\nu_k$, $r_k$ are to be determined later. $G(\eta)$ is a solution of the following differential equation:

where $\delta_k$, $k = 0, 1, 3, 4$, are given in Table 2. In all the cases mentioned above and the subsequent solutions, I used the mathematical software Mathematica to find the values of the constants $\nu_0, \, \nu_1, \, \nu_2, \, r_1, \, r_2$ and $w$. Thus, the analytic solutions to (3.5) using the generalized algebraic method will be presented here with different values of the constants $\delta_k$, $k = 0, 1, 3, 4$.

$(5).$ When $\delta_0 = \dfrac{\delta^2_2}{4 \delta_4},$ $\delta_1 = \delta_3 = 0,$ $\delta_2 < 0$, $\delta_4 > 0$, and $\varepsilon = \pm1,$

$\textbf{Case 7.}$ The solution is given by

Figure 3 shows the time evolution of the analytic solutions. Figure 3$(a)$ shows $\Psi_7$ with $t$ = $0:2:6$. The parameter values are $\delta_2 = -1$, $\delta_4 = 1,$ $\epsilon = -1,$ $\alpha = 0.50,$ $\beta = 0.6,$ $\gamma = -1.5,$ $\lambda = 1.8$ and $x0 = -10.$ Figure 3$(b)$ shows $\Psi_8$ with $t = 0:2:8$. The parameter values are $\delta_2 = 1$, $\delta_4 = -1,$ $\epsilon = -1,$ $\alpha = 0.50,$ $\beta = 0.6,$ $\gamma = -1.5,$ $\lambda = 1.8$ and $x0 = -10.$ Figures 4–6 present the $3D$ time evolution of the analytic solutions $\Psi_2$ (left) and the numerical solutions (right) obtained employing the scheme 5.1 with $t = 5, 15, 25$, $M_x = 1600,$ $N_y = 100,$ $x = 0\to60$ and $y = 0\to1.$

$(6).$ When $\delta_0 = 0,$ $\delta_1 = \delta_3 = 0,$ $\delta_2 > 0$, $\delta_4 < 0$, and $\varepsilon = \pm1,$

$\textbf{Case 8.}$ The solution is given by

$(7).$ When $\delta_0 = 0,$ $\delta_1 = \delta_4 = 0,$ $\delta_3\neq0,$ $\delta_2 > 0$, $\varepsilon = \pm1$

$\textbf{Case 9.}$The solution is given by

$\textbf{Case 10.}$ The solution is given by

4.   Numerical solution

In this section I extract numerical solutions to the resulting ODE system (3.5) using several numerical methods. The purpose of this procedure is to guarantee the accuracy of the analytic solutions. I picked one of the analytic solutions above to be a sample, (3.11). The nonlinear shooting and BVP methods, at $t = 0,$ are used by taking the value of $\psi$ at the right endpoint of the domain $\eta = 0$ with guessing the initial value for $\psi_\eta.$ The new target is to obtain the second boundary condition of $\psi$ at the left endpoint of a particular domain. Once the numerical result is obtained, I compare it with the analytic solution (3.11). The MATLAB solver ODE15s and FSOLVE ^[47] are used to get the numerical solution. The resulting ODE (3.5) is discretized as

for the BVP method and

for the shooting method. Figure 7 presents the comparison between the numerical solutions obtained using the above numerical methods and the analytic solution. Figure 7 shows that the solutions are identical to the analytic solution.

Thus, it is possible to verify the correctness of the analytic solution. I also accept the obtained numerical solution as an initial condition for the numerical scheme in the next section.

5.   Numerical scheme on a fixed mesh

In this section, I use the finite-difference method to obtain the numerical results of system (1.1) over the domain $[a, \, b]\times [c, \, d].$ Here, $a, \; b, \; c$ and $d$ represent the endpoints of the rectangular domain in the $x$ and $y$ directions, respectively, and $T_f$ is a certain time. The domain $[a, \, b]\times [c, \, d]$ is split into $(M_x+1)\times (N_y+1)$ mesh points:

where $\Delta_x$ and $\Delta_y$ are the step-sizes of the $x$ and $y$ domains, respectively. The system (1.1) is converted to an ODE system by discretizing the space derivatives while keeping the time derivative continuous. Completing this yields

where

subject to the boundary conditions:

Equation (5.2) permits us to use fictitious points in estimating the space derivatives at the domain's endpoints. The initial conditions are generated by

where $\alpha, \, \beta, \, \gamma$ and $\lambda$ are user-defined parameters. In all the numerical results shown in this section, the parameter values are fixed as $\alpha = 0.50, \, \beta = 0.6, \, \gamma = -1.50, \, \lambda = 1.80, \, x0 = -45.0$, $y = 0\to1,$ $x = 0\to60$ and $t = 0\to 25.$ The above system is solved by using an ODE solver in FORTRAN called the DDASPK solver ^[48]. This solver used a backward differentiation formula. Since I do not have the initial conditions for the space derivatives, I approximate the Jacobian matrix of the linearized system by using LU-Factorization. The obtained numerical results are acceptable. This can be observed from the Figures 8 and 9.

6.   Stability of the numerical scheme

The von Neumann analysis is used to examine the stability of the scheme (5.1). The von Neumann analysis is occasionally called Fourier analysis and is utilized exclusively when the scheme is linear. Hence, I suppose that the linear version is given by

where $s_0 = \gamma \Phi$, $s_1 = \gamma \Psi$, $s_2 = \lambda \Gamma$, $s_3 = \lambda \Psi$ are constants. Since $\Gamma_y = \Psi_x,$ and $\Psi_y = \Phi_x,$ the first equation of (6.1) is given by

where $\alpha, \, \beta, \, \gamma, \, \lambda, \, s_0, \, s_1, \, s_2, \, s_3, \, l_4$ are constants. I set directly

and also I can have

Substituting (6.3) into (6.2) and doing some operations, I have

Hence,

where

Thus,

The stability condition of the von Neumann analysis is fulfilled. Consequently, from Eq (6.5), the scheme is unconditionally stable.

7.   Error analysis

To examine the accuracy of the numerical scheme (5.1), I study the truncation error utilizing Taylor expansions. Suppose that the error is

where $\Psi(x_m, y_n, t_{k+1})$ and $\Psi_{m, n}^{k+1}$ are the analytic solution and an approximate solution, respectively. Substituting (7.1) into (5.1) gives

where

Hence,

Accordingly, the truncation error of the numerical scheme is

8.   Results and discussion

I have prosperously employed several analytical methods to extract the traveling wave solutions to the two-dimensional Novikov-Veselov system, confirming the solutions with numerical results obtained using the numerical scheme (5.1). The major highlights of the results are shown in Table 3 and Figures 8–10, which allow immediate comparison of the analytic solutions with the numerical results. Through these, I can notice that the solutions are identical to a large extent, and the error approaches zero whenever the value of $\Delta_x, \Delta_y\to 0$. The numerical schemes are unconditionally stable for fixing the parameter values $\alpha = 0.50, \, \beta = 0.6, \, \gamma = -1.50, \, \lambda = 1.80, \, x0 = -45.0$, $y = 0\to1,$ $x = 0\to60$ and $t = 0\to 25.$

Figure 1 presents the time evolution of the analytic solutions $(a)$ $\Psi_1$ and $(b)$ $\Psi_2$ with $t = 0, 10, 20$. The parameter values are $x0 = -20, \ \alpha = 0.50, \ \beta = 0.6, \ \gamma = -1.5,$ and $\lambda = 1.$ Figure 2 presents the wave behavior by changing a certain parameter value and fixing the values of the others. Figure 2$(a, b)$ presents the behavior of $\Psi_1$ when I change the values of (a) $\alpha$ or $\beta$ and (b) $\gamma$ or $\lambda$. In Figure 2$(a)$ it can also be seen that the value of $\alpha$ or $\beta$ affects the direction and amplitude of the wave, such that a negative value always makes the wave negative, its amplitude decreases when $\alpha, \, \beta \to 0$, and its amplitude increases when $\alpha, \, \beta \to \infty.$ In Figure 2$(b)$ the value of $\gamma$ or $\lambda$ affects the direction and amplitude of the wave, such that a negative value always makes the wave negative, and its amplitude decreases when the value of $\gamma$ or $\lambda$ increases. In Figure 2$(c, d)$ present the wave behavior of $\Psi_2.$ Figure 3 shows the time evolution of the analytic solutions. Figure 3$(a)$ shows $\Psi_7$ with $t = 0:2:6$. The parameter values are $\delta_2 = -1$, $\delta_4 = 1,$ $\epsilon = -1,$ $\alpha = 0.50,$ $\beta = 0.6,$ $\gamma = -1.5,$ $\lambda = 1.8$ and $x0 = -10.$ Figure 3$(b)$ shows $\Psi_8$ with $t = 0:2:8$. The parameter values are $\delta_2 = 1$, $\delta_4 = -1,$ $\epsilon = -1,$ $\alpha = 0.50,$ $\beta = 0.6,$ $\gamma = -1.5,$ $\lambda = 1.8$ and $x0 = -10.$ Figures 4–6 present the $3D$ time evolution of the analytic solutions $\Psi_2$ (left) and the numerical solutions (right) obtained employing the scheme 5.1 with $t = 5, 15, 25$, $M_x = 1600,$ $N_y = 100,$ $x = 0\to60$ and $y = 0\to1.$ These figures provide us with an adequate answer that the numerical and analytic solutions are quite identical. Barman et al. ^[42] accepted several traveling wave solutions for (1.1) as hyperbolic functions. The authors employed other parameters to develop new forms for the accepted solution. They proposed that Eq (1.1) describes tidal and tsunami waves, electromagnetic waves in transmission cables and magneto-sound and ion waves in plasma. In comparison, I have found numerous solutions also as hyperbolic functions. Furthermore, I obtained the numerical solutions to enhance the assurance that the solutions presented here are correct and accurate.

9.   Conclusions

I have successfully utilized the generalized algebraic and modified F-expansion methods to acquire the soliton solutions for the two-dimensional Novikov-Veselov system, verifying these solutions with numerical results obtained by employing the numerical scheme (5.1). The major highlights of the results shown in Figures 8–10 and Table 3, which allow immediate comparison of the analytic solutions with the numerical results. Through these, I can notice that the solutions are identical to a large extent, and the error approaches zero whenever the value of $\Delta_x, \Delta_y\to 0$. The numerical schemes are unconditionally stable for fixing the parameter values $\alpha = 0.50, \, \beta = 0.6, \, \gamma = -1.50, \, \lambda = 1.80,$ $x0 = -45.0$, $y = 0\to1,$ $x = 0\to60$ and $t = 0\to 25.$ The Jacobi elliptic functions have effectively deteriorated to hyperbolic functions. The applied numerical schemes have provided reliable numerical solutions when using a small value of $\Delta_x, \Delta_y\to0$.

Ultimately, I can deduce that the methods used are valuable and applicable to extract soliton solutions for other nonlinear evolutionary systems found in chemistry, engineering, physics and other sciences.

Conflict of interest

The author declares that he has no potential conflict of interest in this article.