Influence of numerical discretization on chaos: a quantitative study of Lyapunov dynamics and bifurcation behavior

1 Introduction

Chaos describes the behavior observed in non-linear dynamic systems that is determined and unpredictable [14, 31]. Even though chaotic phenomena exhibit randomness, there is some form of hidden order, which can be explained using mathematics [13, 30]. The most distinctive feature of chaotic systems is their extreme sensitivity to initial conditions, which is popularly known as the butterfly effect, since small changes produce exponentially diverging paths and preclude predictions beyond a certain time span [30]. Most of these systems cannot be solved analytically, hence numerical simulation is vital in understanding their dynamics. Numerical techniques represent continuous equations by breaking them down into discrete steps, although truncation and rounding errors can affect the results dramatically. For instance, using the Euler technique to solve the differential equation y′(t) = cos(t) results in approximating the analytical solution y(t) = sin(t), at discrete intervals. In the case of chaotic systems, errors from numerical discretization can either suppress real chaos or create fake chaos [22].

The difficulty inspires the quest for reliable techniques which would be able to detect chaotic numerical methods which do not alter the dynamics after extended periods of integration. There are known classical fixed-step methods, advanced adaptive algorithms, and particular oscillatory techniques based on trigonometric polynomials. As it is known from literature, any general numerical algorithm may impose stability or instability effects foreign to the original dynamical system [20]. Thus, assessing the reliability of numerical simulations necessitates chaos quantifiers such as phase spaces, periodograms, Lyapunov exponents, and Kolmogorov-Sinai entropy.

There exist numerous uses of chaos theory in the areas of cryptography, biology, climate modeling, economics, astronomy, and neurology [18, 27]. Since there are quite a lot of chaotic dynamical systems exhibiting oscillations, one might apply techniques suggested by Gautschi [24], exploiting frequency behavior. Therefore, the current investigation will attempt to prove the feasibility of the use of such numerical algorithms in preserving the chaotic behavior after extended time periods.

1.1 Theoretical framework

Chaos Theory is one of the fields under non-linear dynamical systems that deals with deterministic systems which show high sensitivity to initial conditions [19]. A typical continuous autonomous dynamical system can be represented by:

dxdt=f(x),x(t0)=x0

For the linearized system

ξ=J´(xn)ξ*

The stability matrices of Numerical schemes are

MFE=I+hj, MBE=(I-hj)-1Numerical stability (Mh)≤1,λ1(h)> 0

where the future course of development rests solely on the current position. Chaotic systems are typified by their boundedness, non-periodicity, strange attractors, and exponentially diverging orbits [26]. The phenomenon of sensitive dependence is measured via the Lyapunov exponent:

λ=limN→∞1N∑n=1N-1ln∣f′(xn)∣

The Lyapunov exponent can be considered positive for a chaotic system, while K-S entropy serves as a sufficient condition for the same when it is positive [23]. Phase space reconstruction allows estimating chaotic invariants from observations of time series [34]

λ1(h)=limN→∞(1Nh)∑n=0N-1(h)∥Dϕh(xn)∥

Where λ1(h)>0 , ∑iλi(h)< 0

Numerical simulations are necessary because there are few analytical solutions to chaotic systems. Time is quantized using numerical methods to discretize continuous-time trajectories. The most basic method is the Euler method:

yn+1=yn+hfn

Whereas, higher-order Runge-Kutta schemes increase both local accuracy and stability [10]. Nevertheless, discretization leads to truncation and round-off errors which either dampen genuine chaos or introduce fake chaotic dynamics. According to the shadowing lemma, even though numerical computations are subject to errors, their solutions are still close to real solutions [24]. In the case of oscillatory chaotic dynamics, Gautschi [35] introduced trigonometric polynomial techniques where the coefficients depend on the oscillation parameter v = ωh:

∑j=0kαjyn+j=h∑j=0kβj(v)fn+j,v=ωhβ0(v)→-12 , β0(v)→32 as ν→0

The proposed approaches will ensure minimum phase error and enhance long-term stability for oscillatory dynamics. The following representative chaotic models, namely, the Lorenz, Rossler, Chen, Chua, Rabinovich–Fabrikant, and Lorenz–Stenflo models, were considered due to their non-linear and oscillatory dynamics.

For the continues system

x=f.(x)

An a Numerical method define as a discrete flow map

xn+1=ϕh(xh)

The method is continuing if

ϕh(xh)=x+hf(x)+Oh(p+1)

1.2 Literature review

A wide range of applications of chaotic models includes but is not limited to neural networks, cryptography, signal processing, optimization, population biology, chemistry, and fluid dynamics [28]. However, simulating the behavior of chaotic models has always been problematic, since previous research focuses more on practical applicability of the numerical approaches than ensuring their correctness for modeling chaotic dynamics. Chaotic models demonstrate sensitivity to initial conditions, therefore, any slight deviations in numerical solutions might develop exponentially over time. Clean numerical simulation increases numerical accuracy through high-order Taylor series and ultra-precise computations; however, supercomputers are required for performing simulations [35].

From simulation investigations, it is found out that discretization errors have substantial effects on chaotic behavior. Truncation and round-off errors might reduce real chaos or induce fake chaotic behavior [12, 29]. Backward Euler methods tend to inhibit chaos when large integration time steps are employed [15, 32], while forward Euler and Runge-Kutta algorithms might lead to spurious chaos under unstable discretization settings [16, 33]. Long-term simulations continue to pose difficulties due to the frequent occurrence of numerical errors that disrupt integration processes.

The theoretical basis was enriched by integrating fundamental works on numerical integration and quantification of chaos, as well as classical references on ODE solvers, Lyapunov exponent computation and long-time numerical behavior of chaotic systems [1].

However, most simulation methods use Euler and Runge-Kutta integration schemes, but their validity regarding chaos preservation remains unproven. Other methods such as Adomian decomposition, differential quadrature, spectral and collocation procedures show good short-time behavior but lack testing for their ability to maintain chaotic invariants over a longer time, e.g., through evaluation of Lyapunov exponents and Kolmogorov Sinai entropy [2]. Another technique, trigonometric polynomials introduced by Awrejcewicz et al. [3], seems to be promising for systems that show oscillations, as they make use of natural frequency characteristics. Systematic comparison with the standard adaptive algorithms over long simulation times is rare. Therefore, there remains an open problem to find.

Although significant progress has been achieved in chaos theory and numerical simulation, there is still a major research gap with regard to the development of computationally effective numerical techniques that can effectively retain chaotic dynamics during long-time simulations without introducing any artificial chaos or chaos dampening. Research to date focuses predominantly on well-established Euler- and Runge-Kutta-type numerical algorithms. In addition, simulation tests have mostly been performed over a short period of time, and their accuracy in terms of rigorous indicators of chaos such as the Lyapunov exponent and Kolmogorov-Sinai entropy has been tested infrequently. Moreover, few comparisons have been made between oscillatory numerical algorithms based on trigonometric polynomials. To fill this research gap, the current paper proposes the formulation of an approach for preserving chaos in non-linear oscillators using trigonometric polynomial-based numerical methods.

Originality of this work is in the quantitative comparison of classical, adaptive and Gautschi-type trigonometric schemes for the long-term preservation of chaos in terms of Lyapunov exponents, K-S entropy, spectral behavior and computational cost.

In the past, the numerical artifacts of chaotic systems have been explored, but most work has been performed using conventional Euler- or Runge-Kutta-type schemes, short duration simulations or single chaos indicators. However, a systematic long-time comparison of fixed-step, adaptive and trigonometric-fitted methods based on common chaos-preservation criteria is still missing. The present work differs from the previous studies in that it investigates both numerical chaos suppression and false chaos generation, and compares the performance of backward Euler, forward Euler, AB–AM, ODEX, RKF45 and Gautschi-type schemes in preserving chaotic invariants for long simulations using phase portraits, periodograms, Lyapunov exponents, K-S entropy and computational cost.

2 Methodology

Chaotic systems have proved their usefulness especially in cryptography, whereby image encryption algorithms utilize the chaotic system's high dependence on initial conditions. For example, the Chen chaotic system, which is numerically simulated using the Runge-Kutta method, is used in the generation of private keys for encryption purposes, where the initial conditions define the ciphering and decryption processes [4]. Nevertheless, the accuracy of numerical simulations poses a problem since most chaotic dynamical systems do not have explicit solutions, and numerical methods may introduce different sources of error such as truncation and rounding errors that may suppress chaos or produce false chaos [5, 17]. Therefore, there is a need for a systematic approach in selecting suitable numerical methods that can maintain chaos throughout long periods of simulation. The methodology, which includes three interrelated stages, includes dynamical systems and numerical methods selection, experimental simulation, and qualitative/quantitative chaos analysis.

2.1 Selection of dynamical systems and numerical methods

At this stage, the chosen dynamical systems, integration techniques, and criteria for assessment are selected. Continuous dynamical systems that are described by ordinary differential equations are chosen as the selected dynamical systems. The selected dynamical systems can be further subdivided into two categories. The first category consists of widely-studied chaotic systems such as the Lorenz system, Rössler attractor, Chen system, and double scroll Chua system. Systems that can be chaotic or non-chaotic based on different parameters belong to the second category and include the Lorenz system (different Rayleigh number), Rabinovich-Fabrikant system, and Lorenz-Stenflo system.

The numerical integration techniques have been divided into three categories depending on their mathematical approach and conformity with the features of dynamical systems. Common integration techniques adopt constant integration step size and consist of B–E technique with two corrections and F–E method. They are adopted due to their simple calculations and scientific evidence that backward Euler technique eliminates chaos whereas forward Euler causes false chaos if an inappropriate step size is used. Advanced integration techniques incorporate varying integration step size and order with error control procedures. The methods include R–K F 45, ODEX technique, and AB-AM technique. They are applied for solving general first-order ODE systems. They are adopted from trusted databases Netlib and ACM [6] to avoid rounding errors. Specific integration techniques use fixed integration step size but rely on trigonometric polynomials. Second-order Gautschi integration technique represented by

yn+2-yn+1=h(β1fn+1+β0fn)

with coefficients

β0(v)=-12(1+112v2+1120v4) andβ1(v)=32(1-14v2+ 1120v4)

The Gautschi approach is created for oscillatory second order problems. However, it may be applied to first order chaotic systems in a component-wise form the Used Cautschi type formula

xn+2-xn+1=h[β1(v)f(n+1)+β0(v)f(n)]fV =hw, moreover when v → 0β0(v)→-12 , β1(v)→32

where v = ωh, is particularly suited to oscillatory problems [7]. The oscillatory behavior of chaotic problems is consistent with the method, which assumes that by harnessing their intrinsic oscillatory characteristics, chaos will be maintained for extended periods compared to general-purpose techniques.

2.2 Simulation experiments

The simulation stage conducts three experiments aligned with the mathematical dynamics of chaotic systems. The first experiment establishes chaos suppression cases by simulating second-group systems under chaotic parameter values using backward Euler with varying integration steps. Based on theoretical bounds h_*c*1 from Owolabi et al. [8] and Khater and Temsah [7] when the chosen step h_*s*1>h_*c*1, numerical results are expected to suppress chaos, evaluated qualitatively through phase space. The second experiment establishes false chaos cases by simulating second-group systems under non-chaotic parameter values using forward Euler, where step sizes h_*s*2>h_*c*2 are predicted to generate spurious chaotic behavior. The third experiment aims to preserve chaos in first-group systems over long time intervals [0, 50, 000] using standard (backward Euler only), sophisticated (AB-AM, ODEX, R-K F 45 with error tolerances 10⁻² and 10⁻⁸), and special (Gautschi) methods. The Gautschi method requires frequency approximation for each system: Lorenz (24.0), Rossler (1.5), Chen (24.0), and Chua (3.0), with integration steps determined via v = ωh simulation Experiments.

The domain frequency was estimate from the partial are by w=arg max p(w_*k*)

λ1(h)>0 , ∑iki(h)<0 , ∑iλhi<0

The simulation stage performs three experiments based on mathematical dynamics of chaotic systems. The first experiment develops chaos suppression scenarios by conducting simulations of second group systems with parameters of chaos using backward Euler integration techniques with different integration steps. According to theoretical constraints h_*c*1 obtained from Owolabi et al. [8], if the selected step h_*s*1 > h_*c*1, numerical simulation is predicted to suppress chaos, analyzed quantitatively in phase space. The second experiment creates false chaos scenarios by performing simulations of second group systems with non-chaotic parameters using forward Euler integration technique, where step sizes h_*s*2 > h_*c*2 can be predicted to cause false chaos. The third experiment seeks to sustain chaos within first group systems within large time ranges [0, 50,000] using standard (backward Euler), advanced (AB-AM, ODEX, R-K F 45 with error tolerances 10⁻² and 10⁻⁸), and specialized (Gautschi) techniques. The Gautschi method requires frequency approximation for each system: Lorenz (24.0), Rossler (1.5), Chen (24.0), and Chua (3.0), with integration steps determined via v = ωh simulation Experiments.

2.3 Evaluation metrics and computational implementation

The evaluation step employs both qualitative and quantitative criteria. In the first two tests, the concept of chaos is measured through phase space. The third test uses simulations that have been completed over the whole time interval based on the following criteria: (1) strange attractor phase spaces; (2) periodogram analysis; (3) function evaluation count as a measure of computation cost; (4) Lyapunov exponent by means of lyap_k and lyap_spec routines from TISEAN package [9]; and (5) Kolmogorov-Sinai entropy using TISEAN d2 routine [10].

Lyapunov exponent quantifies the exponential separation of nearby orbits:

λ=limN→∞1N∑n=1N-1ln∣f′(xn)∣

while the K-S entropy is defined as:

K=limτ→0limδ→0limL→∞1Lτ(KL+1-KL)

A step size convergence test was carried out with

hjho2j, j=0,1,2,……mE(hj)=∥xhj-xh j+1∥2∥xh j+1∥2<ε

when the chaotic indications were steady, i.e.,

⌈λ(hj)−λ(hj+1)⌉< ϵλ,

⌈k(hj)−k(hj+1)⌉<ϵ k

A positive K-S entropy suggests chaotic behavior. Delay coordinate embeddings using Whitney's embedding theorem make possible the estimation of these quantities from one-dimensional time series data, which helps determine if the simulations preserve chaos over long periods.

The Lyapunov exponent and K-S entropy were computed for each of the numerical methods over M independent time blocks. The standard error was calculated as

SE=δM

s is the sample standard deviation. The 95% confidence interval was determined to be

CI95%=x¯ ±1.96 SE

3 Experiments and results

Below, the experiments and results carried out in this research are presented in detail. All experiments were performed using Fortran language.

3.1 First experiment

The objective of this experiment was to identify conditions under which the backward Euler method suppresses chaos. The hypothesis proposed that inappropriate integration steps can eliminate chaotic behavior in numerical simulations.

Chaotic parameter configurations for the second group of dynamical systems were selected from established theoretical studies, together with initial conditions and simulation intervals. Based on the analysis of Varsakelis and Anagnostidis [8], a critical bound h_*c*1 was considered such that chaos is suppressed when the selected integration step satisfies h_*s*1 > h_*c*1, whereas chaotic behavior is preserved for h_*s*1 ≤ h_*c*1. Suitable values of h_*s*1were determined through trial-and-error procedures.

Phase-space analysis was used as the evaluation metric: For the Lorenz system, the parameters were σ = 10, R = 28, and γ = −8/3, with initial conditions (1, 1, 1), simulation interval [0, 50], and bound h_*c*1 = 0.0018. For the Rabinovich–Fabrikant system, α = 1.1, γ = 0.87, initial conditions (−1.0824, 0.1141, 0.2873), interval [0, 50], and h_*c*1 = 0.007 were used. For the Lorenz–Stenflo system, the parameters were σ = 1, r = 26, β = 0.7, and s = 1.5, with initial conditions (1, 1, 1, 1), interval [0, 50], and h_*c*1 = 0.0503.

With the results of the numerical experiments, we conclude that an incorrect choice of the integration step leads to chaos suppression in chaotic systems, so the hypothesis is supported.

The suppression of chaos observed across Figures 1–3 occurs consistently when h_*s*1 > h_*c*1, indicating that backward Euler introduces numerical damping that dominates the non-linear instability. This is not a gradual process, implying that the integration algorithm changes the Jacobian matrix of the system in such a way that the eigenvalues move into the stable area. This means that the chaos suppression effect is caused by the discretization process itself and not inherent in the actual continuous system.

3.2 Second experiment

The purpose of the second experiment was to determine under what conditions the forward Euler algorithm creates false chaos during the numerical simulation of a non-chaotic continuous system. It is assumed that an inadequate integration step will create artificial chaos within a non-chaotic system.

The experimental procedure followed a structure similar to that of the first experiment. Based on previous theoretical studies [11], parameter configurations associated with non-chaotic dynamics in the group 2 test systems were selected together with their initial conditions and simulation intervals.

According to the theoretical analysis, a critical integration threshold h_*c*2 exists such that when the selected integration step satisfies h_*s*2 > h_*c*2, the numerical solution exhibits artificial chaotic behavior. Conversely, when h_*s*2 ≤ h_*c*2, the original non-chaotic dynamics are preserved. The appropriate step sizes were determined through trial-and-error procedures by progressively adjusting the integration step until false chaotic patterns either appeared or disappeared.

Phase-space analysis was used as the principal evaluation metric to distinguish between genuine and numerically induced chaotic behavior. The parameter values, initial conditions, simulation intervals, and critical bounds h_*c*2employed in the experiment are presented below.

Lorenz system: Parameter values are σ = 10, R = 20, γ = −8/3, initial conditions are (1, 1, 1), simulation interval is [0, 50], and the bound is h_*c*2 = 0.0028.

Rabinovich-Fabrikant equations: Parameter values are α = 0.1, γ = 0.9, initial conditions are (1, 1, 1), simulation interval is [0, 100], and the bound is h_*c*2 = 0.009.

Lorenz-Stenflo equations: Parameter values are σ = 1, r = 26, β = 0.7, s = 0.1, initial conditions are (1, 1, 1, 1), simulation interval is [0, 50], and the bound is h_*c*2 = 0.03.

With the results obtained in this experiment Figures 4–6, it is concluded that an inadequate selection of the integration step in the forward Euler method generates false chaos. Consequently, the hypothesis proposed for this experiment is not rejected.

Chaotic attractor formation when h_*s*2 > h_*c*2 signifies a numerical problem associated with artificial divergence. As opposed to genuine chaos, this phenomenon is devoid of consistency and emanates from error magnification. It is hereby evident that forward Euler scheme can exhibit spurious positive Lyapunov characteristics hence generate chaos within supposedly stable systems.

3.3 Third experiment

This third experiment aimed at determining numerical methods that were able to maintain chaos in the test group 1 systems throughout the whole simulation period. It was hypothesized that the special Gautschi scheme was able to maintain chaos within fewer computations as compared to complicated adaptive approaches.

Simulation period was specified as [0, 50,000]. Three categories of numerical methods were analyzed: standard, sophisticated, and special methods.

Among the standard methods, only the backward Euler scheme was considered. Since this method requires a predefined integration step, a trial-and-error strategy was employed, beginning with h = 0.1 and progressively decreasing the step until stable integration over the entire interval was achieved.

For sophisticated methods, AB–AM, ODEX, and Runge–Kutta–Fehlberg 45 were implemented. These methods require relative and absolute tolerances, which were set equal in all experiments using values of 10⁻² and 10⁻⁸. Because these are variable-step methods, an output frequency f_*s* was also defined as an upper limit for the integration step (h ≤ f_*s*). Its value was determined simulation Experiments by progressively reducing f_*s* until the methods completed the full simulation interval.

The special-method category included the Gautschi method. Its integration step depends on the oscillation frequency ω of the chaotic system; therefore, the system frequencies were first estimated numerically. Using these frequencies, the parameter v = ωh was adjusted iteratively, beginning with h = 2⁻¹ = 1/2 and refining it until stable long-term simulations were obtained.

Table 1 reveals a clear relationship between system frequency and the integration step required for stable simulation. High frequency systems, like the Lorenz and Chen systems where ω = 24.0 require much smaller time steps to guarantee numerical stability, whereas lower frequency systems like the Rossler attractor can take larger steps. However, there is an exception to this rule—the Chua system, which suggests that the factors affecting time step requirements are non-linearity and stiffness rather than frequency. Hence, in case of chaotic systems, the time step adjustment based on frequency needs to be accompanied by a stability analysis.

Table 1. Approximate frequencies and corresponding integration steps for chaotic systems using the backward Euler method.

The first set of continuous dynamical systems was selected for this experiment. It should be stressed that all these systems are characterized by chaotic dynamics. Below the parameters of the chosen chaotic systems and initial conditions are provided.

Lorenz system: Parameter values are σ = 10, R = 28, γ = −8/3, and initial conditions are (−15.8, −17.48, 35.64). Rossler attractor: Parameter values are a = b = 0.2, c = 5.7, and initial conditions are (−9, 0, 0). Chen system: Parameter values are a = 35, b = 3, c = 28, with initial conditions (−10, 0, 37). Chua system: Parameter values are α = 10, β = 11.5, m1 = −0.276, m₂ = −3.036, b₂ = 0.8, b3 = 1.4, d₂ = m1b₂, and d3 = d₂ + m₂(b3–b₂), and initial conditions are (0.1, 0, 0).

It is evident from Tables 1, 2 that there is a trade-off between stability and computation in the case at hand. Whereas, AB–AM is able to solve up to the complete time period, ODEX and R-K F 45 break down early on in their calculations, especially when tolerances are tighter. It shows that step-size control does not provide adequate stability in such a system over an extended period.

Table 2. Numerical results of basic, advanced and Gautschi schemes in chaotic problems.

The unified dataset (Table 2) clearly shows that the accuracy of the method depends both on the tolerance of errors and on the form of the used algorithm. Although the method AB–AM maintains the chaotic behavior irrespective of the tolerance level, its stability is provided at the expense of intensive evaluation. In turn, ODEX and R-K F 45 are unstable methods incapable of sustaining the process of calculation for a long time, especially under more rigorous tolerances. The Gautschi method can be considered exceptional because it is a stable and effective approach capable of performing full interval integration with much greater step sizes.

3.3.1 Evaluation and comparison

The metrics used to evaluate the numerical results of this experiment are the following: phase space and periodogram qualitatively; number of evaluations, Lyapunov exponents, and K-S entropy quantitatively.

3.3.2 Phase space

The attractor consistency across Figures 7–10 indicates that most numerical methods preserve the global geometric structure of chaotic dynamics. However, distortions observed in lower-precision AB–AM simulations suggest that insufficient tolerance leads to topological deformation of attractors, which may affect invariant measures and long-term statistical properties.

3.3.3 Periodogram

The broadband spectral distribution observed in Table 3 confirms the presence of deterministic chaos across all systems, characterized by continuous frequency spectra rather than discrete peaks. Nevertheless, the difference that appears in the case of the Chua circuit when employing AB-AM with EEr=10-2 shows a lack of accuracy, which causes the loss of high frequencies. Therefore, it can be stated that accuracy plays a crucial role in obtaining a correct spectrum of the system.

Table 3. Plots for periodograms demonstrating similar behavior to white noise due to the presence of high-frequency peaks, which is typical of chaotic systems; the Chua circuit using AB-AM with Er = 10⁻² up to 500 s.

3.3.4 Number of evaluations

Tables 4, 5 demonstrates that the Gautschi technique leads to a considerable gain in terms of computational expense when compared to backward Euler and high-precision AB-AM integration techniques. In Figure 11 the savings in terms of function evaluation amount to close to two orders of magnitude without compromising the dynamics of the system, as verified by Lyapunov exponent and K-S entropy calculations. These gains arise from the ability of the Gautschi scheme to adapt its frequency to match that of the oscillations of the chaotic trajectory.

Table 4. Number of function evaluations in chaotic systems employing various numerical integration techniques.

Table 5. Lyapunov exponent spectrum estimates showing λ1 > 0 and ∑λ_i < 0 for nearly all systems, confirming chaos, except the Chua system (first 500 s only due to lyap_spec limitations).

^*First 500 s of numerical simulation were taken into account.

In addition to the number of function evaluations, the CPU execution time and the peak memory usage were also measured for each method using the identical hardware configuration. The computational efficiency was assessed using

η=NfTCUP

3.3.5 Lyapunov exponents

Positive maximum Lyapunov exponents indicate the presence of chaos in all cases, suggesting persistent chaos. Comparison among the various approaches reveals that quantifying chaos is significantly influenced by the discretization process. The Gautschi approach offers robust and statistically significant results with reasonable variance, implying successful retention of dynamic instability.

Positive K-S entropies for all systems indicate chaotic behavior. For most scenarios, the Gautschi approach yields larger entropy values, demonstrating better preservation of information production and oscillation dynamics. Thus, the use of frequency adaptive integration allows better capturing the unpredictability inherent in chaos dynamics.

To verify the methods generated Lyapunov exponents and attractor structures were compared with published benchmark results of the Lorenz, Rossler, Chen and Chua systems. Agreement in the sign of, positivity of K-S entropy and preservation of the known attractor geometry were utilized as proof that the numerical method preserves the expected chaotic dynamics. Additional validations could be performed on other benchmark systems such as the Duffing oscillator, Hindmarsh–Rose model and hyperchaotic Lorenz-type systems.

Thus, in all scenarios considered here, the methods discussed maintain chaos over long intervals of simulations, with the exception of backward Euler that fails to do so with improper step sizes. Overall, Tables 5, 6 the Gautschi approach seems to provide the best compromise between computation and dynamics preservation. It preserves chaotic patterns with minimal number of evaluations, along with positive Lyapunov exponents and K-S entropies.

Table 6. K-S entropy estimates for chaotic systems using different numerical methods with standard errors, 95% confidence intervals, and positivity-confirmed chaotic behavior (first 500 s only).

^*First 500 s of numerical simulation were taken into account.

It was shown in the first two experiments that backward Euler could suppress chaos, while forward Euler could create chaos when inappropriate integration steps were used. As shown in the third experiment, chaos could be maintained for long intervals using the approaches mentioned above—backward Euler, AB–AM with relative error 10⁻⁸, and Gautschi techniques over extended periods of time.

On the other hand, Gautschi demonstrated greater efficiency and was able to generate phase spaces, periodograms, and statistical parameters at least comparable, if not superior, to those generated by the other algorithms.

The third important feature of Gautschi's algorithm is that, from the second iteration onward, only one calculation per step is needed. As such, the algorithm's computational complexity would be similar to Euler methods while maintaining greater dynamics accuracy. Moreover, when the oscillation frequency, or its approximation, is known, the selection of the appropriate step for integration becomes relatively easy.

Adaptive approaches control the local truncation error, however in chaotic systems small local mistakes are exponentially magnified, as per

∥e(t)∥≅C tal eλti

Thus, local error control does not ensure long-time preservation of chaotic invariants. On the other hand, the Gautschi technique uses frequency-dependent coefficients through v = wh which decreases the phase error and enhances the long-time stability for oscillatory chaotic trajectories.

4 Discussion

The results show that numerical discretization significantly affects the qualitative and quantitative dynamics of non-linear chaotic systems. Instead of being neutral approximating algorithms, numerical integration can hinder, support, or enhance chaos based on their stability properties and discretization criteria. This finding expands previous research on numerical chaos with a systematic study comparing several methods in their ability to reproduce chaotic behavior.

Numerical dissipation was evident in the results of the backward Euler scheme. Chaotic orbits tended to damp and converge to fake fixed points once the discretization step passed a certain limit. This finding was confirmed with the analysis of the Lorenz system and Rabinovich-Fabrikant oscillator, proving that implicit schemes could mislead the actual dynamics due to over-stabilization. The forward Euler integration, on the other hand, created false chaos for large integration steps. Exaggerated rounding errors resulted in chaotic and even attractor-type orbits despite the absence of inherent chaotic dynamics in those systems.

Comparison of advanced numerical methods showed large discrepancies regarding the preservation of chaos over long periods of time. The Adams Bashforth-Adams Moulton (AB-AM) approach [11] was able to preserve the chaotic regime when applied with very strict tolerances; however, it demanded extremely expensive computation. The ODEX and Runge-Kutta-Fehlberg integrators often failed at completing long-term runs, which implies that adaptive local error control cannot ensure global stability.

The Gautschi algorithm turned out to be the most efficient among all considered methods. In case of the Lorenz model, it took more than 10 million iterations to simulate the chaotic dynamics using the AB-AM scheme with tolerance 10⁻⁸; however, the Gautschi scheme managed to complete the process with 300 thousand iterations, thus reducing computation cost almost by 97 percent. Similar results were achieved for the Chua model.

Lyapunov exponents indicated that all three schemes maintained positive divergence, although the divergence magnitude was highly dependent on the scheme used. With regards to the Lorenz system, the AB–AM method yielded results close to the expected value, while Gautschi was found to have underpredicted the divergence rate because of phase approximation. On the other hand, with the Chua attractor, Gautschi showed significantly higher Lyapunov exponents, indicating the tendency of the numerical scheme to exaggerate chaos for highly non-linear systems.

Moreover, Kolmogorov-Sinai entropy showed how numerical methods can affect information generation rates. Gautschi often resulted in higher entropy values than traditional methods, implying more efficient representation of oscillation behavior, albeit with larger entropy values potentially indicative of numerical exaggeration. Comparative analysis across the two attractors showed that Lorenz dynamics remained fairly consistent with discretization, while the Chua attractor had high sensitivity due to multiple scroll formation.

In conclusion, the study highlights the influence of numerical integration on chaotic invariants and not the other way around. The developed assessment tool, combining considerations for efficiency, Lyapunov divergence, and entropy production, serves as a reliable foundation for measuring numerical accuracy in non-linear dynamic systems [25]. The discovery of the Gautschi method as a fast and accurate solver algorithm will have great significance in the field of applied mathematics and beyond.

5 Conclusions and future work

This research examines the issue of chaos persistence when numerically simulating oscillating non-linear systems over extended periods. Chaotic systems are very sensitive to small changes in their initial states, non-periodic, and usually have no analytical solution, which makes the process of numerical approximation hard. To explore this question, a three-step approach involving selection, simulation, and evaluation stages was applied. Several representative chaotic systems of dynamical processes and numerical integrators like Gautschi method, Euler, Runge–Kutta Fehlberg, ODEX, and Adams Bashforth–Moulton were compared. The emphasis on the Gautschi method is determined by the nature of the oscillating dynamical system. Chaos persistence was quantitatively measured based on theoretically substantiated criteria, such as Kolmogorov-Sinai entropy, Lyapunov exponents, and phase space analysis (Table 7). It has been found out that the choice of the method has a critical impact on chaos persistence. The research proposes a consistent foundation for the selection of numerical solver for chaotic systems and can be extended for non-chaotic ones, exponential fitting, and cryptographic purposes.

Table 7. Maximum Lyapunov exponent values for chaotic systems employing different numerical schemes with their standard errors and 95% confidence levels.

^*First 500 s of numerical simulation were taken into account.

The primary quantitative conclusion is that the Gautschi scheme retained positive Lyapunov exponents and positive K-S entropy throughout the whole 50,000 s simulation period, despite needing much fewer function evaluations than AB–AM and backward Euler. For example, in the Lorenz system Gautschi used 300,000 evaluations compared to 10,944,766 for high precision AB–AM and 25,000,000 for backward Euler. This shows its practical significance for long-time simulation of oscillatory chaotic systems. The key constraint is the reliance on reliable frequency estimate, and future work should explore multi-frequency fitting, adaptive Gautschi versions, and more chaotic benchmarks.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Sagar R. Khirsariya, Marwadi University, India

Reviewed by: Rakesh Kumar Meena, MIT Vishwaprayag University, India

Anil Chavada, Maharaja Sayajirao University of Baroda, India