Solving a system of nonlinear difference equations with bilinear dynamics

1.   Introduction

Difference equations have long been integral to the analysis of discrete dynamical systems, providing essential tools across fields such as economics, biology, engineering, and time series analysis (see, ^[1,2]). Since the 18th century, significant advancements have been made in deriving closed-form solutions for specific classes of these systems, establishing a robust foundation for exploring their solvability—a topic that remains highly relevant in contemporary research (see, ^[3,4]). Classical studies have extensively examined the solvability of difference equations, focusing on exact solutions for specific, well-defined systems (see, ^[5,6]). Additional contributions have further expanded this area, offering new perspectives on foundational principles (see, ^[6,7]). Recent advances, however, have extended the field's scope, tackling more complex systems that challenge traditional analytical techniques. For example, research has examined the global behavior of rational difference equations (see, ^[8,9]) and nonlinear systems (see, ^[10,11]). Moreover, studies on systems involving co-balancing numbers and periodicity have provided novel insights into their dynamic properties (see, ^[12]). Analytical methodologies have been developed to investigate higher-order nonlinear systems of difference equations, focusing on their solutions, stability, and numerical simulations (see, ^[13]). Furthermore, innovative approaches have been employed to analyze the global dynamics of systems with exponential-form difference equations, providing deeper insights into their behavior (see, ^[14]).

The study of solvability in difference equations often involves the derivation of closed-form solutions, although these solutions can sometimes become exceedingly complex. In such cases, qualitative approaches—like the examination of invariants and boundedness properties—provide an alternative means to gain insights into system behavior, especially when explicit solutions are either too intricate or infeasible (see, for example, Schinas ^[15], 1997). Despite these challenges, having general solution formulas for new classes of difference equations is invaluable for enhancing our understanding of their underlying dynamics. Often, systems with nonlinear characteristics can be transformed into simpler, solvable forms, allowing them to adopt solvability properties from well-known classes of equations. A notable example of these transformation techniques in the study of difference equations is seen in Stević's research. In his groundbreaking work ^[16],

he introduced innovative mathematical methods that effectively simplified complex two-dimensional systems of difference equations, enabling the derivation of explicit solutions. This achievement has significantly advanced our understanding of solvability in such systems. In a subsequent study ^[17], Stević further demonstrated the power of these transformations by successfully converting the following system,

into a solvable framework. Additionally, in another significant investigation ^[18], he applied a similar transformation approach to the system:

demonstrating how mathematical transformations can provide deeper insights into the system's dynamics. These examples underscore the importance of transformation methodologies in exploring complex behaviors within difference equations. Such techniques not only facilitate the solvability of these systems but also deepen our understanding of the intricate patterns that characterize their behavior.

In addressing the challenges posed by complex dynamical systems, particularly those with high-dimensional or nonlinear behaviors, numerical methods provide essential tools for obtaining approximate solutions when analytic approaches become infeasible. Notably, the second-order backward differentiation formula alternating direction implicit (BDF2 ADI) Galerkin finite element method has been effectively applied to solve three-dimensional evolutionary equations with nonlocal terms. This approach combines spatial discretization via Galerkin finite elements with an ADI scheme, ensuring stability and convergence for the system under specified regularity conditions (Yang et al. ^[19], 2022). Additionally, the super convergence properties of the orthogonal Gauss collocation method (OGCM) offer robust solutions for two-dimensional fourth-order sub-diffusion equations, demonstrating optimal error bounds that maintain robustness as parameters approach critical values (Yang and Zhang ^[20], 2024). These advanced numerical techniques underscore the potential of finite element and finite volume approaches in the analysis and simulation of difference equations, expanding the toolkit available for examining complex system dynamics.

With the rise of increasingly complex systems across a range of scientific fields, the demand for robust models that accurately capture both linear and nonlinear behaviors has grown significantly. Among various classes of difference equations, three-dimensional systems present a particularly rich area for exploration due to their ability to represent real-world phenomena driven by the interaction of multiple variables. This paper focuses on a specific three-dimensional system of difference equations defined as follows:

where $\alpha, \beta, \gamma, \delta, \varepsilon, \lambda, \tau, \sigma\in \mathbb{R},$ $\alpha^{2}+\beta^{2}\neq0, \gamma^{2}+\delta^{2}\neq 0, \varepsilon^{2}+\lambda^{2}\neq0, \tau^{2}+\sigma^{2}\neq0,$ and $P_{0}, Q_{0}, P_{-1}, Q_{-1}, R_{0}\in\mathbb{R}.$

Modeling dynamic systems through difference equations has long proven to be an effective approach for analyzing complex behaviors across disciplines, from biology and economics to physics and engineering. In particular, recursive difference equations offer a powerful framework for examining the evolution of interconnected variables over time. Such systems capture both short-term fluctuations and long-term trends that arise from the interplay between different factors, making them especially suitable for applications in biological rhythms, social dynamics, and environmental interactions.

This paper examines a three-dimensional nonlinear system of difference equations designed to model fundamental behavioral dynamics within animal populations. The system relies on three primary variables: $P_{m},$ $Q_{m}$, and $R_{m}$, which represent individual behavior, social interactions, and environmental pressures, respectively. Through recursive relations, the model captures how these factors evolve in response to internal and external stimuli, offering a comprehensive perspective on the system's dynamics. By transforming the system into a two-dimensional bilinear form, the analysis becomes more tractable while maintaining the intricate interactions among the variables.

The main objective of this study is to investigate the solvability of the system and to derive closed-form solutions under various parameter conditions. These theoretical findings are enriched by extensive numerical simulations that visualize the dynamic trajectories of the variables $P$, $Q$, and $R$ over time. The simulations not only confirm the analytical results but also reveal the complex interplay between individual actions, social interactions, and environmental factors.

Furthermore, this model has practical applications in understanding animal behavior under changing environmental conditions. By simulating key behavioral metrics—such as movement patterns, interaction frequencies, and responses to environmental stressors—the model provides a structured framework for analyzing how animal populations adapt to fluctuations in their environment. This approach is especially valuable for studying ecosystem stability and adaptability, as it offers predictive insights into how various factors influence behavioral dynamics over time.

The significance of this research lies in its integration of recursive difference equations with dynamic simulations, creating a flexible and cohesive framework applicable to a wide range of scenarios, from ecological systems to engineered environments. This study not only contributes to the theoretical foundations of nonlinear systems but also opens pathways for future research into more complex models and multifaceted interactions.

The structure of this paper is organized as follows: Section 2 presents a detailed analysis of the specific system under investigation, including the derivation of closed-form solutions. Section 3 focuses on the simulations, providing a dynamic analysis of animal behavior and system stability using recursive models. Finally, Section 4 concludes the paper by summarizing the findings and opening avenues for future applications.

2.   Solving the equation system in (1.a)

In this section, we delve into a systematic methodology to solve the difference equations presented in (1.a), considering specific initial conditions. The sequence $\left(P_{m}, Q_{m}, R_{m}\right)$ represents the solution to the system. A critical observation here is that the system can collapse or lose its definition if any initial values are set to zero. To preserve the solution's consistency and avoid singularities, it is imperative that the condition $R_{m+1}P_{m}Q_{m}\neq0$ holds for $m\geq-1$. This condition ensures the continued existence of a viable solution. From this point forward, we assume that $R_{0}P_{-k}Q_{-k}\neq0$ for $k\in\left\{ 0, 1\right\}$, which is crucial to maintaining the integrity of the system over time. With these precautions, we can explore the system's dynamics and uncover potential behaviors that emerge as the system evolves.

To facilitate further analysis, we introduce a change of variables to simplify the relationships among the system's variables. This transformation provides a fresh perspective on the system and is defined as follows:

This transformation restructures the system into a more manageable form, setting the stage for extracting key properties of the model. After applying this change, the system (1.a) transforms into the following form:

This new formulation reveals a bilinear structure, which simplifies the analysis and provides deeper insights into the interactions between $\widehat{P}_{m}$ and $\widehat{R}_{m}$. With this bilinear form, we can explore scenarios where the system demonstrates stability, periodicity, or even chaotic behavior, depending on the initial conditions and parameter values chosen.

The two-dimensional nonlinear system of difference equations, now reformulated as (2.b), was originally introduced by Stević and Tollu ^[18]. In their work, they demonstrated the solvability of this system by offering a comprehensive approach for finding its general solution. Their analysis highlighted that specific parameter constraints are necessary for the system's solvability. For example, when $\gamma = 0$, the first equation of (2.b) can be substituted into the second equation, resulting in a simplified structure characterized by interlacing indices—a configuration derived from a bilinear equation. This special case becomes solvable due to the reduction in complexity. Similarly, a corresponding simplification occurs when $\tau = 0$, producing a dual case to the one where $\gamma = 0$. Given these observations, we will operate under the assumption that neither $\gamma$ nor $\tau$ is zero, ensuring that the system retains its nonlinear complexity. Moreover, to avoid ill-defined solutions, we impose the additional conditions that $\gamma\widehat{R}_{m}+\delta\neq0$ and $\tau\widehat{P}_{m}+\sigma\neq0$ for all $m$, ensuring the terms in the equations remain valid throughout iterations.

With these assumptions established, we present the following lemma, based on the foundational work of Stević and Tollu ^[18], which forms the cornerstone for further investigation of the system.

Lemma 2.1. Let $\alpha, \beta, \gamma, \delta, \varepsilon, \lambda, \tau, \sigma\in\mathbb{R},$ where the following condition: $\gamma\left(\alpha^{2}+\beta^{2}\right) \tau\left(\varepsilon^{2}+\lambda^{2}\right) \neq0$ holds$.$ Additionally, let $\rho_{1}$ and $\rho_{2}$ represent the roots of the quadratic polynomial:

Under these conditions, the system of difference equations given in (2.b) admits a closed-form solution. The general solution to this system can be described in two distinct cases, depending on whether the roots $\rho_{1}$ and $\rho_{2}$ coincide or differ.

Case 1. $\rho_{1} = \rho_{2}.$ In the case of a repeated root $\rho_{1}$, the general solution to the system is given by:

Case 2. $\rho_{1}\neq\rho_{2}.$ When the roots $\rho _{1}$ and $\rho_{2}$ are distinct, the general solution reflects a more complex interplay between the two roots:

In Lemma 2.1, the general solution to the system of difference equations is derived based on whether the roots $\rho_{1}$ and $\rho_{2}$ of the characteristic polynomial $\Lambda\left(\rho\right)$ are identical or distinct. This approach allows us to categorize the solution into two main cases, each reflecting different dynamics of the system:

Case 1. Repeated root ($\rho_{1} = \rho_{2}$). When the roots are identical, the general solution exhibits a linear growth behavior. Here, each term in the solution sequence depends on a combination of initial conditions and the repeated root $\rho_{1}$. This form typically involves terms that increase linearly with the index $m$, reflecting a progressive change in the system without oscillatory components. The structure of the solution in this case is particularly straightforward, as it follows a predictable pattern influenced primarily by the initial values and $\rho_{1}$.

Case 2. Distinct roots ($\rho_{1}\neq\rho_{2}$). When the roots are distinct, the solution becomes more complex, showing a mix of exponential or oscillatory behavior depending on the magnitudes and signs of $\rho_{1}$ and $\rho_{2}$. The terms of the solution sequence will include expressions involving powers of $\rho_{1}$ and $\rho_{2}$, leading to dynamic behavior that can alternate, grow, or decay over time. This case captures the interaction between the two roots, revealing richer dynamics compared to the repeated root case.

The general solution in both cases depends heavily on the values of the parameters and initial conditions. By categorizing the solution this way, we can analyze the system's behavior under various scenarios, providing insights into its stability, oscillatory nature, and long-term behavior.

Remark 2.1. In Lemma 2.1, the condition $\gamma\left(\alpha^{2}+\beta ^{2}\right) \tau\left(\varepsilon^{2}+\lambda^{2}\right) \neq0$ is essential for ensuring the existence and nature of the roots of the quadratic polynomial $\Lambda\left(\rho\right)$. This condition serves several critical purposes. First, it guarantees that the parameters produce real or complex roots, which are necessary to express the system's closed-form solution. Without this condition, the polynomial might yield roots that are not suitable for constructing a solution, thereby complicating or even preventing the solution process. Second, the condition enables us to analyze the system's long-term behavior effectively. When the roots are distinct and real, the system might exhibit oscillatory or divergent behavior, while repeated roots lead to a different, often linear, growth behavior in the solution. This distinction in behavior is foundational for understanding the system dynamics. Lastly, the condition helps avoid singularities in the solution. Certain parameter values might otherwise cause undefined expressions or singularities in the solution, especially where denominators approach zero. Therefore, the condition ensures that the system remains mathematically tractable and that the solution is well-defined for all relevant parameter values.

Remark 2.2. In certain special cases, the system can be further simplified when specific relationships hold among the parameters. For instance, when $\beta \gamma = \alpha\delta$ or $\lambda\tau = \varepsilon\sigma$, the system (2.b) reduces significantly. In the case of $\beta\gamma = \alpha\delta$, the sequences $\left(\widehat{P}_{m}\right)$ and $\left(\widehat{R}_{m}\right)$ stabilize to constant values, expressed as:

Similarly, when $\lambda\tau = \varepsilon\sigma$, the sequences $\left(\widehat{P}_{m}\right)$ and $\left(\widehat{R}_{m}\right)$ are reduced to constants governed by a different set of parameters:

These simplified cases are invaluable when certain parameter combinations lead to degenerate behavior in the system, allowing us to deduce constant sequences for $\left(\widehat{P}_{m}\right)$ and $\left(\widehat{R}_{m}\right)$, which facilitates further analysis.

Remark 2.3. In our analysis of the quadratic polynomial $\Lambda\left(\rho\right)$, the condition

guarantees the generality of the solution by ensuring that the roots, $\rho_{1}$ and $\rho_{2}$ remain distinct. If this condition were not met, the roots would coincide, leading to a single repeated solution structure that limits flexibility in analyzing and solving the system's dynamics under different initial conditions. Thus, maintaining distinct roots allows for a broader exploration of the system's behavior across various cases and parameter variations, which is especially valuable in dynamic analysis and closed-form solution derivations, as shown in the specific solutions given in Lemma 2.1 for cases of both distinct and repeated roots.

According to Lemma 2.1, we derive the closed-form expressions for the solutions of system (1.a). Building on these solutions, we can now formulate the following theorem.

Theorem 2.1. Let $\alpha, \beta, \gamma, \delta, \varepsilon, \lambda, \tau, \sigma\in\mathbb{R},$ where the condition $\gamma\left(\alpha^{2}+\beta^{2}\right) \tau\left(\varepsilon^{2}+\lambda^{2}\right) \neq0$ holds$.$ Under these parameters, the system of difference equations presented in (2.b) admits a general closed-form solution. The solution to this system is expressed in two distinct cases, based on whether the roots $\rho_{1}$ and $\rho_{2}$ are identical or distinct.

Case 1. When the assumption in (2.c) is not satisfied, the general solution becomes:

for $s\in\left\{ 0, 1\right\},$ and $m\geq0,$ where $\widehat{P}_{0} = \dfrac{P_{-1}}{P_{0}Q_{0}}$ and $\widehat{R}_{0} = \dfrac{R_{0}}{P_{-1}}.$

Case 2. Under the assumption specified in (2.c), the general solution to the system is:

for $s\in\left\{ 0, 1\right\},$ and $m\geq0.$

Proof. To prove this theorem, we begin by transforming the system (1.a) into a more tractable form. This is achieved by introducing a set of variable substitutions, as follows:

From the second equation of the system (2.d), we find:

Additionally, by separating terms based on parity, we obtain:

Substituting this expression for $Q_{m}$ into the first equation of the system (2.d), we derive an explicit form for $P_{m}$:

Further, by dividing terms into even and odd cases, we have:

Finally, using the expression for $P_{m}$ from (2.f), and applying the change of variables from (2.a), we derive that: $R_{m} = P_{m-1}\widehat{R}_{m}$ for $m\geq0.$ By breaking down further, we get: $R_{2m} = P_{2m-1}\widehat{R}_{2m}$ and $R_{2m+1} = P_{2m}\widehat{R}_{2m+1}$ for $m\geq0$. □

Example 2.1. In this example, we consider a system defined by three sequences: $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$. The evolution of these sequences follows the recursive relations given by system (1.a), where $\alpha = 0.56$, $\beta = -0.5$, $\gamma = 0.9$, $\delta = 1.2$, $\varepsilon = 1.5$, $\lambda = -0.86$, $\tau = 0.78$, and $\sigma = 0.46$. The initial values are set as $P_{-1} = 0.8$, $P_{0} = 0.5$, $Q_{-1} = 0.4$, $Q_{0} = 0.6$, and $R_{0} = -0.3$. To visualize the behavior of these sequences over multiple iterations, we generate plots for system (1.a), as shown in Figure 1.

In this example, the recursive sequences $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$ exhibit intriguing dynamics as illustrated in Figure 1. The sequence $\left(P_{m}\right)$ displays mild oscillations that suggest a gradual stabilization toward a steady state, reflecting the system's tendency to find equilibrium over time. Conversely, $\left(Q_{m}\right)$ demonstrates pronounced oscillatory behavior, characterized by substantial fluctuations that indicate a heightened sensitivity to the initial conditions and the recursive relationships governing the system. Meanwhile, $\left(R_{m}\right)$ presents a more restrained oscillation pattern, highlighting its balanced interaction with the other sequences. The combined analysis of the even and odd indexed sequences for $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$ reveals significant differences in their oscillatory characteristics, emphasizing the complex interplay within the system and its ongoing evolution toward convergence.

In Figure 1, we initially observe the dynamics of the sequences $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$ produced by the recursive model, capturing fundamental patterns like oscillations and tendencies toward stabilization. However, while Figure 1 offers valuable insights into the behavior of the recursive system, it is essential to validate its accuracy and reliability by comparing it with established numerical methods. This transition to Figures 2 and 3 allows for a comparative analysis that helps assess the efficacy of our model against the Newton and Runge-Kutta methods, both known for their precision in solving differential and difference equations. Figures 2 and 3, therefore, provide critical validation of our approach by illustrating how closely the outputs from our recursive method align with those of the Newton and Runge-Kutta methods. This comparative analysis enables a thorough evaluation of the model's robustness, ensuring that it captures the essential dynamics and reliably approximates the results produced by conventional numerical techniques. By highlighting both the strengths and any limitations of the recursive model in Figures 2 and 3, we ensure that the proposed method is suitable for complex system analysis and understand where it might require adjustments for improved precision.

In Figure 2, we compare the dynamics of the sequences $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$ generated by the recursive approach against those produced by the Newton method. This comparison shows a close alignment between both methods across all three sequences, indicating that the recursive model effectively approximates the Newton method's results. For $\left(P_{m}\right)$, both methods capture similar oscillatory patterns, though the Newton method shows slightly smoother transitions. The $\left(Q_{m}\right)$ sequence exhibits notable oscillations, with both methods mirroring these fluctuations closely, albeit with minor differences in peak values. Similarly, $\left(R_{m}\right)$ shows consistency between the two approaches, demonstrating that the recursive model effectively replicates the stability and oscillations of the Newton method. Overall, Figure 2 suggests that the recursive method provides a reasonable approximation to the Newton method, making it a viable computational alternative.

In Figure 3, we compare the trajectories of the sequences $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$ under the recursive model with those generated by the Runge-Kutta method. This comparison shows that both methods produce similar oscillatory patterns, with each sequence displaying trends toward stabilization over time. For $\left(P_{m}\right)$, the recursive model closely follows the Runge-Kutta trajectory, especially in later iterations, suggesting a convergence toward steady-state behavior. For $\left(Q_{m}\right)$, the recursive approach captures the Runge-Kutta method's oscillations and sensitivity with high fidelity, although occasional minor deviations in amplitude are observed. Similarly, $\left(R_{m}\right)$ shows a strong alignment across both methods, highlighting the recursive model's reliability in maintaining accuracy over multiple iterations. Figure 3 thus supports the robustness of the recursive model as it approximates the Runge-Kutta method, especially in capturing cyclical behavior and stability trends, underscoring its practical use in complex system analysis.

Example 2.2. In this example, we consider a system defined by three sequences: $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$. The evolution of these sequences follows the recursive relations given by system (1.a), where $\alpha = 0.90$, $\beta = -1.10$, $\gamma = 0.80$, $\delta = 1.40$, $\varepsilon = 1.70$, $\lambda = 0.90$, $\tau = -0.70$, and $\sigma = 0.85$. The initial values are set as $P_{-1} = 1.20$, $P_{0} = -1.10$, $Q_{-1} = 0.60$, $Q_{0} = 0.95$, and $R_{0} = 0.90$. To visualize the behavior of these sequences over multiple iterations, we generate plots for system (1.a), as shown in Figure 2.

In this example, the system defined by the sequences $\left(P_{m}\right)$, $\left(Q_{m}\right)$, and $\left(R_{m}\right)$ demonstrates pronounced oscillatory dynamics, as depicted in Figure 4. The sequence $\left(P_{m}\right)$ exhibits heightened oscillations compared to precedent example, indicating increased sensitivity to changes in the recursive structure and interactions with other sequences. $\left(Q_{m}\right)$ stands out with considerable amplitude oscillations, underscoring its dynamic and unstable nature, which appears to be influenced significantly by both initial conditions and the recursive relationships. Conversely, $R_{m}$ shows consistent oscillatory behavior with values that tend to stabilize, albeit with less dramatic fluctuations than $Q_{m}$. The analysis of the even and odd indexed sequences for each variable further highlights the system's complexities, where distinct patterns emerge, revealing a system that remains highly dynamic and far from equilibrium. Overall, the sequences showcase a robust interdependence that drives the ongoing evolution of the system.

3.   Dynamic animal behavior modeling: A comprehensive analysis of system simulations

In this section, we examine animal behavior dynamics through a simulation model based on a system of recursive difference equations. These equations capture the temporal evolution of essential behavioral metrics, including social interactions, environmental pressures, distances traveled, and overall behavioral patterns. By simulating real-world scenarios, this model provides a structured framework to analyze both short-term fluctuations and long-term trends. The recursive nature of the system enables us to track how each variable—individual behavior, group dynamics, or environmental influences—evolves in response to external stimuli and internal interactions. This modeling approach is particularly useful for assessing system stability, adaptability, and the impact of external pressures, ultimately offering predictive insights into future behaviors under different conditions.

The system is driven by three main components represented by the variables $P,$ $Q$, and $R$. The variable $P$ reflects individual behavior or position, evolving based on both past behavior and current environmental conditions. $Q$ represents social interactions, capturing the influence of group dynamics, such as cooperation or competition, on individual actions. $R$ enotes environmental pressures, such as resource availability, which dynamically interact with both animal behavior and social interactions. The evolution of these variables is governed by recursive difference equations, enabling us to model the interdependencies between individual actions, social dynamics, and environmental factors. Through this approach, we can explore how animal behavior adapts over time, gaining insights into the underlying patterns that influence system stability and adaptability.

This graph illustrates the temporal evolution of behavioral metrics $P$, $Q$, and $R$ over 50 time iterations. The behaviors demonstrate oscillatory dynamics, indicating fluctuations in the animals' responses. The gradual changes in amplitude suggest the presence of both short-term reactions and long-term adjustments to their environment. The distinct patterns seen in $P$, $Q$, and $R$ highlight their individual responses to external stimuli. This transition to more stable states might indicate a shift toward equilibrium, signifying how the system balances out over time. This figure is foundational for understanding the overall system dynamics, capturing key oscillations and trends in animal behavior. It sets the stage for the transition to Figure 6, which explores social interactions in greater depth. Social interactions are essential to the system, as they both influence and are influenced by the behavioral dynamics observed in Figure 5, providing a more comprehensive view of the interconnected factors driving animal behavior.

Figure 6 illustrates the frequency of social interactions among the animals throughout the simulation. The graph shows intermittent spikes in social activity, followed by periods of lower interaction. These fluctuations may result from changes in environmental pressures or individual behavioral shifts. Notably, an upward trend emerges, indicating that social interactions become more frequent as the system evolves, suggesting a growing intensity in group dynamics over time. This upward trend aligns with the behavioral changes observed in Figure 5, underscoring the influence of social interactions on the model's dynamics. The next step, shown in Figure 7, is to examine the spatial behavior of the animals, specifically focusing on the distance traveled between them over time. Moving from social interactions to spatial movement highlights how these interactions often drive physical movement, whether animals come together or drift apart.

This graph illustrates the spatial behavior of the animals, focusing on the distance between them over time. The fluctuating distances indicate dynamic movement patterns, with periods of closeness followed by separation. A closer look reveals that as social interactions increase (Figure 6), the animals tend to either approach each other or maintain a varying distance. This suggests a correlation between social interactions and spatial distribution, where the need for either interaction or avoidance influences the distances observed. Understanding these spatial dynamics is crucial for interpreting movement behaviors in real-world ecological systems. Building on this, Figure 8 introduces environmental pressure, a critical external factor influencing both social interactions and spatial movements.

In this figure, environmental pressure is plotted as a function of time. The graph reveals subtle oscillations, indicating that the environment exerts a consistent yet variable influence on the animals' behaviors. This ongoing pressure shapes the adaptation of $P$, $Q$, and $R$ over time, particularly in their responses to both internal and external stimuli. The stability of the environmental pressure curve reinforces the notion that while the environment plays a significant role, the primary drivers of change within the system may stem more from social interactions and spatial dynamics, as illustrated in Figures 6 and 7. To integrate these complex relationships, Figure 9 presents the $3D$ trajectory of $P$ vs. $Q$ and $P$ vs. $R$, showcasing how these variables interact within multidimensional space. The transition to this $3D$ visualization is motivated by the need to observe how individual behavior, social dynamics, and environmental influences intertwine over time.

Figure 9 presents a three-dimensional visualization of the trajectories for $P$ vs. $Q$ and $P$ vs. $R$ over time, illustrating the complex interactions between these variables. Both trajectories show circular or spiral patterns, indicating cyclical behaviors within the system. In the $P-Q$ plot, we observe a gradual convergence toward a more stable region, suggesting that the system could eventually achieve a dynamic equilibrium. Similarly, the $P-R$ trajectory exhibits a comparable stabilization trend. These $3D$ plots provide a clear depiction of the interconnections among the system's components, highlighting how $P$, $Q$, and $R$ influence one another across multiple dimensions. The transition from this view of trajectories to the next figure, which examines the spatial distribution of animal positions, is essential for gaining a comprehensive understanding of not only their movement patterns but also the stability of animals within their environment.

Figure 10 presents the distribution of animals' positions over time, represented by the variables $P$, $Q$, and $R$. The plot shows that these positions span a wide range, indicating diverse movement patterns among the animals. Notably, $R$ exhibits higher variability compared to $P$ and $Q$, which display more moderate fluctuations. This variation suggests that some animals are more mobile, continuing to shift locations significantly, while others remain relatively stable, occupying certain areas over time. This distribution offers a clear picture of the animals' spatial dynamics, shedding light on how they disperse and organize within the environment modeled. Some groups appear to gravitate toward specific regions, while others exhibit ongoing movement. This spatial behavior is critical for understanding how individual and group positioning evolves over time. Moving to Figure 11, we examine the relationship between individual behavior and group dynamics, building on the insights from Figure 10's trajectories and distributions to further explore the interconnectedness of these behaviors.

Figure 11 illustrates the relationship between the variables $P$ and $Q$, highlighting their correlation across time. The linear trend of the trajectory suggests a direct, proportional relationship between $P$ and $Q$, indicating that fluctuations in $P$ are closely mirrored by corresponding changes in $Q$. The clustering near the origin reflects that, as the system approaches stability, the values of $P$ and $Q$ converge toward specific points. This plot supports the observations from Figure 9, reinforcing the interdependence between these variables. Next, the heatmap provides a more comprehensive view of activity within the system, visually distinguishing periods of high and low engagement. This broad perspective on dynamic interactions complements the more focused analyses in earlier figures, offering an immediate visual summary of the system's overall behavior and variability.

In Figure 12, the heatmap provides a detailed visualization of the dynamics of variables $P$, $Q$, and $R$ across time. The color gradients represent the intensity of each variable's values throughout the simulation, where darker shades signify higher values and lighter shades indicate lower values. This color-coded approach allows for an immediate understanding of periods of intense activity or stability in the animals' behavior. The heatmap highlights correlations and variations among the variables, illustrating how $P$, $Q$, and $R$ fluctuate in response to each other and to environmental conditions. This visualization captures the interactions and cyclical patterns that may not be as evident in other plots, offering a comprehensive view of how each behavior metric evolves over time. Transitioning to a figure depicting stability regions for parameters $\alpha$ and $\beta$ is essential, as it links the behavioral trends and interactions observed here to the parameters that shape the system, providing insights into the stability and adaptability of the modeled environment.

Figure 13 depicts the stability regions for parameters $\alpha$ and $\beta$, offering insights into the conditions that determine whether the system remains stable or transitions into instability. The contour plot visually separates stable regions from unstable ones, with lighter shades indicating areas of stability and darker shades representing zones of instability. This visual analysis is key to understanding the influence of $\alpha$ and $\beta$ on the system's dynamics, as it reveals how specific parameter adjustments can drive the system toward or away from stability. Such insights are fundamental for identifying parameter ranges that maintain balanced behaviors and for predicting the effects of external influences on the overall system.

4.   Conclusions

In this paper, we have conducted an in-depth analysis of the dynamic properties of a three-dimensional nonlinear system of difference equations, systematically reduced to a two-dimensional bilinear form for a more focused investigation. Through recursive relations in the sequences $P_{m}$, $Q_{m}$, and $R_{m}$, representing individual behaviors, social interactions, and environmental pressures, we captured complex dynamics and provided general closed-form solutions that elucidate the stability and adaptability of systems under varying parameter conditions. Extensive numerical simulations complemented our theoretical findings, offering robust insights into the evolution of these variables over time.

Future work could enhance our framework by considering the preservation of essential physical properties, such as maintaining discrete maximum principles $(DMP)$ and ensuring bounded solutions across distorted meshes. Notably, recent advancements in nonlinear finite volume ($NFV$) schemes, such as the work by Yang and Zhang (^[21], $2024$) on $DMP$-preserving approaches for two-dimensional sub-diffusion equations on distorted meshes, open promising avenues for extending this system's applicability. Such methods prevent spurious oscillations and uphold physical constraints on key quantities, which is crucial for modeling realistic behaviors in complex environments. Incorporating these strategies may improve the physical fidelity of models based on recursive difference equations and broaden their relevance to fields where physical constraints are paramount, paving the way for future studies to explore the influence of intricate interactions and environmental variability on system dynamics.

Author contributions

Hashem Althagafi: Methodology, validation, formal analysis, investigation, visualization, writing review, and editing; Ahmed Ghezal: Software, validation, resources, data curation, writing-original draft preparation, supervision, project administration. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

We would like to express our gratitude to the Editor-in-Chief of the journal, the Assistant Editors, and the anonymous referees for their constructive comments, valuable suggestions, and remarks, which were instrumental in improving the final version of the paper. Additionally, we extend our thanks to Prof. Dr. Mouhamed Ghezal and our colleagues, Prof. Dr. Karam Ghezal and Prof. Dr. Rahaf Ghezal, for their important encouragement and support.

Conflict of interest

All authors declare no conflict of interest in this paper.