Analytical Study of Nonlinear Systems of Higher-Order Difference Equations: Solutions, Stability, and Numerical Simulations

: This paper aims to derive analytical expressions for solutions of fractional bidimensional systems of difference equations with higher-order terms under specific parametric conditions. Ad-ditionally, formulations of solutions for one-dimensional equations derived from these systems are explored. Furthermore, rigorous proof is provided for the local stability of the unique positive equilibrium point of the proposed systems. The theoretical findings are validated through numerical examples using MATLAB, facilitating graphical illustrations of the results.


Introduction
Difference equations and systems constitute a foundational element of mathematical modeling, exerting significant influence across a spectrum of scientific domains and practical applications.Over the past decade, this field has witnessed a notable evolution, characterized by a surge in scholarly interest and activity.This trend underscores an increasingly recognized importance and utility of difference equations and systems in tackling intricate phenomena and facilitating informed decision-making processes.
In recent years, there has been a growing emphasis on exploring nonlinear systems of difference equations, driven by a quest for analytical solutions, deeper insights into dynamic behaviors, and applications spanning diverse disciplines including biology, physics, probability theory, environmental science, and engineering.These systems serve as discrete analogs to differential equations, offering a robust framework for modeling time-series data such as economic indicators (e.g., Gross Domestic Product, inflation rates, exchange rates) inherently measured at discrete intervals (see, [1,2]).Existing research has probed various facets of nonlinear dynamics, encompassing local dynamics, topological classifications, bifurcation analysis, and chaos control.For example, Khan et al. [3] investigated these dynamics within the context of a discrete-time COVID-19 epidemic model, while other studies [4,5] have explored similar phenomena across different domains.However, while several methods exist for solving linear difference equations, the landscape of nonlinear systems remains largely uncharted.Despite recent attempts to simplify complex nonlinear systems into linear forms, there persists a significant gap in analytical approaches for addressing systems of difference equations, posing a challenge for researchers striving to deepen their understanding of their behavior and properties (see, [6,7]).
The heightened interest in nonlinear systems of difference equations underscores their potential to capture intricate dynamics and phenomena that may elude representation by linear models.Nonlinear systems often exhibit rich behavior, including bifurcations, chaos, and complex attractors, rendering them valuable tools for comprehending the dynamics of real-world phenomena.A particularly active area of research centers on the pursuit of closed-form solutions to nonlinear systems of difference equations, to enhance our understanding of dynamic systems and refine predictions with greater precision.(see, [8][9][10][11][12][13]).Many researchers have made significant contributions to unraveling the behavior of solved difference equations and systems, shedding light on their dynamics and stability properties.For instance, Elsayed and Alshabi [14] delved into the solution forms and stability properties of second-order systems, while Al-Basyouni and Elsayed [15] provided formulas for solutions to systems of rational difference equations of various orders, demonstrating periodicity in certain cases.Okumuş and Soykan [16] extended this exploration to two-dimensional systems associated with Tribonacci numbers.Our paper extends beyond these realms by delving into the analysis of higher-order systems and considering more diverse and complex scenarios.By addressing specific systems of difference equations and investigating their existence, form, and stability properties, our paper offers a comprehensive examination of nonlinear dynamics in discrete time.Furthermore, we can obtain numbers that satisfy some conditions of interest, for example, Richard and Raoul numbers.Moreover, our paper builds upon the contributions of Ghezal and Zemmouri [12] and Abo-Zeid et al. [6,7] by synthesizing and advancing existing knowledge in the field of nonlinear difference equations and systems.By incorporating recent advancements in theoretical analysis and numerical validation techniques, our study pushes the boundaries of understanding in this domain and opens up new avenues for exploration.While Abo-Zeid [17] provides explicit solutions for specific difference equations, our paper explores the fundamental properties of systems of difference equations representing different analytical challenges.Our work provides a deeper analysis of the behavior of bidimensional systems, while Abo-Zeid's [17] work revolves around studying the details of explicit solutions and the general behavior of specific models.The study by Simşek et al. [18] focuses on solving a specific rational difference equation, with a particular focus on providing explicit solutions and analyzing the behavior of the solutions.In contrast, our paper delves into the examination of a broader range of rational difference equations, exploring various forms and complexities.
In comparing our paper to existing research, we identify opportunities to build upon current knowledge and extend the frontier of understanding in this field.Through rigorous theoretical analysis and numerical validation, our paper aims to fill critical gaps in the literature by addressing specific systems of difference equations and investigating their existence, form, and stability properties.By elucidating the dynamics of nonlinear systems in discrete time, our work seeks to advance the understanding of complex mathematical models and their applications across various scientific and engineering domains.Specifically, in this paper, we focus on exploring the existence and form of solutions for the systems of Difference Equations ( 1)-( 4).
The remainder of this paper is structured as follows: Section 1 offers a detailed review of the literature, highlighting key developments and identifying gaps that motivate the current study.In Section 2, we present the theoretical framework and methodologies employed in our analysis, laying the groundwork for subsequent discussions.Section 3 delves into the empirical findings, presenting numerical simulations and graphical illustrations that corroborate the theoretical insights.Finally, Section 4 synthesizes the key findings, discusses their implications, and outlines avenues for future research.

Main Results
In this paper, we delve into the examination of the subsequent systems of rational difference equations: endowed with initial conditions denoted as σ −2 , σ −1 , σ 0 , ρ −2 , ρ −1 , and ρ 0 , ensuring that they are real and non-zero.
∀n ≥ 0, , q ≥ 0, endowed with initial conditions denoted as σ −j , and ρ −j , j ∈ {0, 1, . . . ,3q + 2}, ensuring that they are real and non-zero.To ensure the robustness and validity of the systems under consideration, the denominators must remain non-zero to prevent any division by zero errors.Hence, a critical condition governing these denominators is that the expression within them should never equate to zero for all n ≥ 0. By imposing this condition, we guarantee the stability and solvability of the system, thereby facilitating a comprehensive and accurate analysis of its behavior and properties.The physical background of the systems of rational difference equations described in Equations ( 1)-(4) lies in their representation of dynamic processes influenced by feedback mechanisms.These systems model scenarios where quantities or variables interact with each other over discrete time intervals, and their behavior depends on their past values and certain parameters.System (1) represents a feedback loop where two variables, σ and ρ, influence each other's evolution over time.The evolution of each variable depends on its past values and the other variable's value at the previous time step.This system could describe phenomena in various fields such as population dynamics, economic systems, or chemical reactions in which two entities interact in a feedback loop.System (2) extends the concept of (1) by introducing a delay parameter q, indicating that the influence between σ and ρ occurs over multiple time steps.This delay could represent phenomena in which there is a time lag between the cause and effect in the system, such as in delayed feedback control systems or systems with transport delays.System (3) introduces nonlinearities in the feedback loop by incorporating quadratic terms in the evolution equations of σ and ρ.This nonlinearity can capture more complex dynamics compared to linear feedback systems.It may be relevant in systems in which interactions between variables exhibit nonlinear behavior, such as in biological systems or nonlinear control systems.System (4) extends the nonlinear feedback model of (3) by introducing a delay parameter q, similar to (2).This combination of nonlinearities and delay further enriches the dynamics of the system, allowing for the modeling of more intricate feedback mechanisms with time delays.These systems of rational difference equations provide mathematical frameworks for understanding the behavior of dynamic systems subject to feedback interactions and time delays, with applications ranging from natural and social sciences to engineering and control theory.
In the context of this study, we thoroughly examine the well-defined solutions of the previously mentioned systems, elucidating their forms and behaviors.This analytical scrutiny is carried out through the introduction of specific substitutions to the phase variables: Remark 1.The term 'well-defined solutions' for Systems (1)-( 4) refers to solutions derived from the sets F = {σ −2 , σ −1 , σ 0 , ρ −2 , ρ −1 , and ρ 0 }, selected outside the respective singularity sets of the systems.A set of F that produces a solution {(σ n , ρ n ), n ≥ 1} for (1), wherein at least one denominator becomes zero at some least index n, leads to an undefined value for σ n+1 and/or ρ n+1 .This particular set is commonly referred to as the 'forbidden set,' as discussed further in [19].

Formulation of Solutions for System (1)
To derive the closed-form solutions for the system represented by (1), we employ the substitution defined in (5).Following some transformations, we arrive at the subsequent system of linear difference equations: By manipulating System (6), we derive the subsequent system ∀n ≥ 0, which transforms into a system of two independent linear difference equations ∀n ≥ 0, This transformation is achieved through the following change of variables: ∀n, ε n = ω n + π n and τ n = ω n − π n .Therefore, the exact closed-form of the general solution for System (8) (resp.System ( 6)) is derived, as detailed in the following Lemma 1 (resp.Lemma 2).
Lemma 1.Let {(ε n , τ n ), n ≥ −3} be the solution to System (8), which comprises homogeneous linear difference equations with constant coefficients along with initial conditions ε and τ 0 ∈ R. For all n ≥ 0, the solutions are given by ∀n, Proof.In solving the first (respectively, second) linear difference equation of System (8), we typically utilize the characteristic polynomial: The roots of this equation are (respectively, where α 1 and α 2 are a pair of complex conjugate roots, and α 3 is a real root satisfying the following conditions: Hence, the closed form of the general solution for the first (respectively, second) linear difference equation of system (8) is given by: where α 1 , α 2 and α 3 are real roots satisfying: , and τ 0 are initial conditions such that: ), and we have where the sequences f j , j = 1, 2, 3, 4 , g j , j = 1, 2, 3, 4 are solutions of the latter systems.
After some calculations, we obtain: The lemma is proven.
Through the above discussion and leveraging Lemma 2, we readily derive the closed form of the general solution for System (1), as presented in the following theorem: Theorem 1.Consider {(σ n , ρ n ), n ≥ 0} as a solution to System (1).Then, for all n, , where, for j = 1, 2, 3, Proof.Utilizing the change of variables (5), the relationships are established as follows: By Lemma 2, the remaining steps are straightforward and, therefore, omitted.
Theorem 3. The positive equilibrium point Ξ is locally asymptotically stable.
Proof.For the proof, we define functions M 1 , M 2 : (0, +∞) 2(3q+2) → (0, +∞) as follows: where c ′ 0:q = c 0 , c 1 , . . ., c q .Linearizing System (2) around the equilibrium point Ξ is a common technique to facilitate its study.To achieve this, we introduce vector c In light of these symbols, we acquire the ensuing representation: where , with O k×l is the matrix with zero entries and I m×m represents the identity matrix.Following some initial computations, the characteristic polynomial of Θ q can be expressed as follows: By utilizing MATLAB R2021b, it is determined that all solutions to Q Θq (λ) = 0, q ≥ 0 reside within the unit disc |λ| < 1.By Rouche's Theorem, it can be concluded that the positive equilibrium point Ξ is locally asymptotically stable.

Formulation of Solutions for System (3)
To derive the closed-form solutions for the system represented by (3), we employ the substitution defined in (5).Following some transformations, we arrive at the subsequent system of linear difference equations: By manipulating System (10), we derive the subsequent system ∀n ≥ 0, which transforms into a system of two independent linear difference equations Therefore, the exact closed-form of the general solution for System (12) (resp.System (10)) is derived, as detailed in the following Lemma 3 (resp.Lemma 4).
After some calculations, we obtain: The lemma is proven.
Through the aforementioned discussion and leveraging Lemma 4, we readily derive the closed form of the general solution for System (3), as presented in the following theorem: Theorem 4. Consider {(σ n , ρ n ), n ≥ 0} as a solution to System (3).Then, for all n, , where, for j = Proof.Utilizing the change of variables (5), the relationships are established as follows: By Lemma 4, the remaining steps are straightforward and therefore omitted.

Remark 3.
In System (4), Ξ constitutes unique positive equilibrium.Now, we introduce a crucial result related to the local stability of System (4).
Theorem 6.The positive equilibrium point Ξ is locally asymptotically stable.
Proof.For the proof, we define functions M 1 , M 2 : (0, +∞) 2(3q+2) → (0, +∞) as follows: Linearizing System (4) around the equilibrium point Ξ is a common technique to facilitate its study.To achieve this, we acquire the ensuing representation: where Following some initial computations, the characteristic polynomial of Θ q can be expressed as follows: 4 , By utilizing MATLAB R2021b, it is determined that all solutions to Q Θ q (λ) = 0, q ≥ 0 reside within the unit disc |λ| < 1.By Rouche's Theorem, it can be concluded that the positive equilibrium point Ξ is locally asymptotically stable.

Conclusions
This paper presents a comprehensive investigation into fractional bidimensional systems of difference equations with higher-order terms, aiming to derive analytical expressions for their solutions under specific parametric conditions.Additionally, formulations of solutions for one-dimensional equations derived from these systems are explored.Rigorous proofs are provided for the local stability of the unique positive equilibrium point of the proposed systems.The theoretical framework established in this study is supported by extensive stability analyses and numerical simulations, which offer valuable insights into the behavior of nonlinear systems of difference equations.Notably, investigating positive equilibrium points and establishing local asymptotic stability for specific systems demonstrate the depth of our understanding.Furthermore, numerical examples using MATLAB are employed to validate the theoretical findings.Graphical illustrations of the results showcase the practical implications of the theoretical insights, reaffirming the stability of the studied systems.However, challenges persist, particularly regarding systems dynamics in unexplored scenarios.This highlights the need for continued research and exploration in this area.This paper lays a robust foundation for future investigations into more complex models, particularly those addressing dynamic systems with non-constant or periodic coefficients.Moreover, it offers avenues for broader applications and potential extensions to various scientific and engineering disciplines.

Figure 1 .Example 2 .
Figure 1.Dynamics exhibited by a specific solution of System (1).Example 2. Figure 2 provides a visual representation of the sustained dynamics within System

Figure 2 .
Figure 2. Dynamics exhibited by a specific solution of System (2).

Figure 4 .
Figure 4. Dynamics exhibited by a specific solution of System (4).