The intricate and often unpredictable nature of chaotic systems, coupled with their prevalence across a wide array of real-world phenomena, has long captivated the attention of researchers, driving extensive and in-depth investigations into their underlying dynamics [1–5]. Recent advancements in the study of chaotic systems have led to the development of increasingly sophisticated analysis techniques, grounded in a robust statistical description framework of the natural measure [6], which weighs the long-term asymptotic behavior on chaotic attractors [7]. The evaluation of averages of various variables, which elucidates macroscopic properties of the system, demands an accurate characterization of the natural measure.

The natural measure is supported on the possibly fractal structure of an attractor, and hence the explicit analytical expression is generally elusive. Unstable periodic orbits (UPOs, also called cycles), instead, may serve as a ‘skeleton' of chaotic dynamics, encapsulating a wealth of dynamical information [8–10]. Cycle expansions, by systematically organizing the set of UPOs in a hierarchical manner, bridge the gap between the global natural measure and local dynamical properties. This technique enables a layered approximation of the natural measure and emerges as a powerful theoretical framework for a quantitative analysis of chaotic systems, when combined with symbolic dynamics [11] and cycle search algorithms [12–17]. Importantly, observable averages can be accurately calculated using only a small number of cycles [18].

In strongly hyperbolic chaotic systems, averages of observables converge exponentially with increasing cycle length due to the shadowing property [7, 19]. However, ubiquitous non-hyperbolicity in physical systems often leads to slow convergence in cycle expansions, limiting their applicability [20]. This usually induces singularities in the natural measure, where the stability eigenvalues of UPOs scale algebraically with topological length (in case of intermittency) near the folding points or their images. Such singularities may not be well captured by a small number of cycles and thus severely degrade the convergence of cycle expansions.

To address this trouble, various techniques have been designed but only make limited improvements, since the singularities are still not accurately quantified. For example, one may [21, 22] derive certain spectral determinant from the distribution of evolution time in hyperbolic regions of intermittent systems, but the accuracy of this approach is poor. Stability ordering scheme [23] may serve as an alternative, converging faster than traditional topological length truncation in non-hyperbolic systems. However, this method requires the stabilities of a large number of long UPOs, and improper stability truncation could lead to huge errors. Infinite symbolic alphabets [24] could be established and used to estimate the fundamental terms of cycle expansion through tail resummation, analytic continuation, or integral transformation. However, this method fails to yield an analytical asymptotic estimate for curvature corrections, leading to substantial numerical errors and limiting its value to purely theoretical interest. The conjugate scheme [25] accelerates the convergence by performing cycle expansion in the conjugate hyperbolic system, while only a very small number of systems have suitable conjugate systems. The induced map [26] reconstructs the natural measures for the singular and non-singular parts using analytic approximation and cycle expansion, respectively, thereby improving computational accuracy and convergence. However, the scheme is relatively computationally intensive and applicable only to one-dimensional intermittent maps. Consequently, based on cycle expansions, it is of paramount importance to develop a program that can be more extensively applied to dynamical systems with a singular natural measure.

The critical issue here is: how to maintain the rapid convergence of the cycle expansion while constructing a computational framework applicable to systems with singular natural measure caused by non-uniform hyperbolicity? In this paper, a novel layered approximation algorithm is proposed to account for the above considerations. It not only addresses the slow convergence of cycle expansions in singular systems but also provides a new perspective for understanding how singularities affect the overall statistical properties, which is of great significance for studying complex non-equilibrium systems.

This paper is organized as follows. Section 2 explores the physical picture, theoretical foundations, and detailed algorithmic aspects of the layered approximation scheme. To verify its effectiveness and generality, section 3 applies it to several typical chaotic systems featuring singular natural measures: those with intermittency, with multiple singularities, or with infinitely many singularities. Finally, section 4 discusses possible improvements and potential applications of the algorithm.

Cycle expansions with finite truncations constitute a hierarchical linear approximation of the dynamics [27]. Local information from densely distributed UPOs approximates the natural measure—a process possibly compromised in non-hyperbolic systems. Slow convergence arises from the failure of the approximation with finite number of cycles to capture singularities in the measure. To address this, we propose a layered decomposition of the natural measure into: (i) singular layers: governed by the distribution of singular points, directly responsible for slow convergence. (ii) Hyperbolic layers: the residual smooth component, being approximated by the natural measure $\tilde{\rho }(x)$ of selected hyperbolic system $\tilde{f}(x)$. $\tilde{f}(x)$ should be topologically conjugated, not necessarily smoothly. The conjugacy is constructed to preserve the symbolic dynamics and the phase-space structure (e.g. the ordering of critical points and their images). Any hyperbolic system whose symbolic dynamics and kneading sequence [7] match those of the original system is acceptable. This decomposition enables rapid convergence in cycle expansions by isolating singularities from hyperbolic dynamics. However, poorly chosen conjugates that alter the kneading sequence may lead to convergence failure.

To validate the hierarchical algorithm and showcase its core principles with maximal clarity and minimal computational cost, we examine a set of prototypical unimodal maps whose natural measures carry distinct singularities. Although one-dimensional, these maps capture concrete physical and biological processes: the logistic-like map governs population dynamics in ecology [28], whereas the intermittent map reproduces the alternation of long laminar intervals and sparse turbulent bursts observed in fluid flows [29]. The Feray-like map has a singularity at x = 0 [as shown in figures 1(a) and (b)]which presents the intermittency as f_1s(x) → x + cx^1+s,  x → 0⁺,  s ∈ (0, 1). The choice of hyperbolic system is non-unique and a practical criterion is to select the simplest piecewise-linear map that (i) preserves the topological entropy and kneading invariants of the original system and (ii) produces a hyperbolic natural measure whose overall shape resembles the singular natural measure with its singular spike removed. The $\tilde{\rho }(x)$ may be derived from the tent map: which has a uniform measure. The logistic map f_2s(x) = 4x(1 − x) has two singularities at x = 0 and x = 1 [as shown in figures 1(c) and (d)], while $\tilde{\rho }(x)$ can also be derived from the tent map (2). The map ${f}_{3s}(x)=\,\rm{sin}\,(\frac{\pi }{a}-\frac{\pi }{a}x),a=1.315644...$ extremely compresses the phase space in the neighborhood of the critical point $x=1-\frac{a}{2}$, making a singularity at x = 1, while $\tilde{\rho }(x)$ may be derived from a linear map: where ${x}_{a}=\,\rm{sin}\,\frac{\pi }{a}=0.684355...$ and x_b = 0.538 770… The image of x = 0 is a fixed point x_a for both maps f_3s(x) and ${\tilde{f}}_{3s}(x)$, which induces three singularities x = 1, 0, x_a in the natural measure, as shown in figures 1(e) and (f). The map ${f}_{4s}(x)=\,\rm{sin}\,(\frac{\pi }{a}-\frac{\pi }{a}x),\,a=1.102634...$ extremely compresses the phase space at $x=1-\frac{a}{2}$ and induces four singularities in the natural measure at x = 1, 0, x_a, x_b, with ${f}_{4s}(1-\frac{a}{2})=1$, f_4s(1) = 0, f_4s(0) = x_a = f_4s(x_b) and f_4s(x_a) = x_b. Clearly, x = x_a and x = x_b constitutes a period-2 cycle $\overline{{x}_{a}{x}_{b}}$. The $\tilde{\rho }(x)$ may be derived from a linear map: where x_a = 0.288 273… , x_b = 0.897 365… , x_d = 0.336 848… , as shown in figures 1(g) and (h).

For any one-dimensional system with a singular natural measure, we may select a hyperbolic map preserving the phase space structure. Singularities arise from the extreme compression of phase space or intermittency, their locations can therefore be inferred directly from the dynamics f(x). The asymptotic form of singularities is analytically derivable from the evolution equation of the natural measure: where x_i = f⁻¹(x) is the ith preimage of x and ${f}^{{\prime} }({x}_{i})$ is the stability factor at x_i. The associated singularity exponent is obtained by matching the orders of x on both sides of the equation. By combining the singular layers with a hyperbolic layer—expanded as a correction function C(x) that acts on the measure $\tilde{\rho }(x)$ of a chosen hyperbolic system whose cycle expansion converges exponentially—we obtain a high-fidelity approximation of the original measure. This enables efficient computation of observable averages well beyond the reach of conventional cycle expansion.

We demonstrate our framework using a prototype system with a singular natural measure (see figure 2), featuring two singularities at x_a and x_b. We can decompose the measure into four singular layers ${\{{\rho }_{i}(x)\}}_{i=1}^{4}$ and one hyperbolic layer ρ₀(x). The asymptotic behavior of each singular layer admits explicit analytical approximation: where the position and the singularity exponent s_i for each singular layer are analytically determined from f(x). For the hyperbolic layer, we write ${\rho }_{0}(x)=C(x)\tilde{\rho }(x)$, where $\tilde{\rho }(x)$ is the natural measure of a conjugate hyperbolic system. Although C(x) is initially unknown, it admits a spectral expansion in terms of orthogonal polynomials. The full natural measure then decomposes as: where the correction function is expanded as $C(x)\,={\sum }_{l}^{n}{c}_{l}{P}_{l}(\varphi (x))$. Away from the corresponding singularity, the singular layer ρ_i(x) is negligibly small and readily cancelled by C(x). {P_l(x)} are the lth order Legendre orthogonal polynomials [30] and ${\{{c}_{l}\}}_{l=0}^{n}$ are the Legendre coefficients. φ(x) = 2x − 1 is a shift-and-scale function and maps phase space interval [0, 1] onto the canonical orthogonal interval [−1, 1] required by Legendre polynomials. The singularity exponents ${\{{s}_{i}\}}_{i=1}^{4}$ can be analytically approximated from the measures in the neighborhoods of the pre-images. The coefficients ${\{{c}_{{s}_{i}}\}}_{i=1}^{4}$ (singular layers) and {c_l} (hyperbolic layer) may be computed numerically. Leveraging the exponential convergence in hyperbolic systems, this transforms the natural measure determination to the evaluation of these coefficients with several exact constraints imposed by: (i) the measure normalization: ${\int }_{{ \mathcal M }}\rho (x){\rm{d}}x=1$. (ii) The dynamical invariance of the natural measure: ⟨a⟩ = ⟨a∘f⟩ for any observable a(x) and the natural measure ρ(x) of the dynamics x_n+1 = f(x_n). This is derived from the conjugacy of the Koopman operator ${ \mathcal K }$ [31] and the Perron–Fröbenius operator ${{ \mathcal L }}_{{\rm{P}}{\rm{F}}}$ [7] under ergodic evolution:

This time-invariance principle, which is intuitive from ergodicity and rigorously established via operator conjugacy, provides critical cues for high-precision measure approximation. Crucially, each observable a_k(x) generates an independent equation: sufficiently many of which is able to determine all $\{{c}_{{s}_{i}}\}$ and {c_l}. The convergence rate of our scheme depends on the smoothness of the chosen observables {a_k(x)}. Smooth observables (e.g. orthogonal polynomials) typically yield faster convergence, while observables with singularities may require more expansion terms The target observable average is then: where the modified observable ${\tilde{a}}_{l}(x)={P}_{l}(\varphi (x))a(x)$ admits efficient computation of $\langle {\tilde{a}}_{l}\rangle $ via rapidly converging cycle expansions [7, 18] (the second term) in the selected hyperbolic system $\tilde{f}(x)$. The second term sums over the combination of all prime cycles and all pseudo-cycles [7]) $\tilde{p}$ of $\tilde{f}(x)$ up to the truncation length ${t}_{{\rm{\max }}}$. Correspondingly, ${A}_{\tilde{p},l}={\sum }_{{x}_{i}\in \tilde{p}}{\tilde{a}}_{l}({x}_{i})$ and ${T}_{\tilde{p}}$, respectively, denote the cumulant observable and the period associated with $\tilde{p}$. The truncation length t_max denotes the maximum topological length of prime cycles or pseudo-cycles retained in the cycle expansion; equivalently, all orbits whose period ${T}_{p}\leqslant {t}_{{\rm{\max }}}$ are included. This cut-off is independent of the expansion order n used for the correction function C(x).

We test the effectiveness of the layered scheme in several example systems and compare it with existing methods. In the numerical experiments, the convergence rates are assessed through a dual-validation mechanism. The averages obtained with the expansion order n + 1 are used as the benchmark values when estimating the logarithmic errors with the expansion order not greater than n. Any observable a(x) can be substituted into our framework and we focused on ⟨x⟩ for clarity and consistency across examples. Moreover, the averages were also computed independently via long-time Monte-Carlo integration [32] to ensure convergence to the correct value. The Monte Carlo validation employs a single trajectory of 10⁹ steps, after discarding 10⁷ transient steps, starting from a randomly chosen typical point x₀ in phase space for each system. Expectation values are computed as

Figures 3(a) and (b) demonstrate accelerated convergence via our layered scheme for the Feray-like map (1). The natural measure at x = 0 (ρ(x) ∼ c_sx^−s as x → 0⁺, $s=\frac{1}{2}$) decomposes into a singular layer and a hyperbolic layer where the tent map (2) is selected as the hyperbolic replacement. Equation (10) can thus be rewritten as: where $\tilde{\rho }(x)$ is the natural measure of the tent map. Here, there is one singular coefficient ${c}_{{s}_{1}}$ and n + 1 hyperbolic coefficients ${\{{c}_{l}\}}_{l=0}^{n}$, totaling n + 2 unknowns. These are determined by applying equation (9) with n + 2 independent observables, e.g. {1, x, x², …, xⁿ⁺¹}. Observable averages converge rapidly with the expansion order n until a critical threshold $\hat{n}$. Beyond $\hat{n}$, the errors increase rather than decrease with the increase of n. Crucially, $\hat{n}$ increases with the truncation length t_max, enabling precision control through parameter selection. Here, it is easy to check from figure 3 that the threshold $\hat{n}$ corresponding to each t_max, is $\hat{n}\approx {t}_{\,\rm{max}\,}-3$. In fact, we may provide a reliable heuristic for parameter selection that $\hat{n}\propto \frac{{t}_{\,\rm{max}\,}}{1.5}$ for the Feray-like map. The relationship is not universal but follows a predictable pattern based on the nonlinearity of f(x). The scaling law $\hat{n}\propto \frac{{t}_{\,\rm{max}\,}}{m}$ depends on the map's local degree m near critical points: m = 1.5 for Feray-like maps, m = 2 for quadratic maps, m = 3 for cubic maps, etc. This exponent m reflects how rapidly the map stretches phase space near singularities. For observables, smoother polynomial bases yield better convergence than non-analytic functions. The relationship is therefore map-dependent but observable-class universal within the same smoothness category. The theoretical explanation for the existence of the threshold $\hat{n}$ is rooted in the term ⟨a∘f⟩ in equation (9), which is employed to determine the expansion coefficients ${\{{c}_{l}\}}_{l=0}^{n}$. The cycle expansion of ⟨a∘f⟩ can only resolve polynomial observables up to a certain order for the dynamics f(x). Given a specific truncation length ${\hat{t}}_{\,\rm{max}\,}$, the most wavy cycle of this length contains $[\frac{{t}_{\,\rm{max}\,}}{1.5}]$ peaks ([… ] denotes the floor function.), which may capture all the wavy profiles of polynomials of order t_max. For this reason, the cycle expansion of ⟨a∘f⟩ can only be effective for a polynomial observable a(x) of order $[\frac{{\hat{t}}_{\,\rm{max}\,}}{1.5}]$ if f(x) is the Feray-like map (order $[\frac{{\hat{t}}_{\,\rm{max}\,}}{2}]$ for quadratic maps, order $[\frac{{\hat{t}}_{\,\rm{max}\,}}{3}]$ for cubic maps and so on). We observed a larger $\hat{n}$ probably because of the existence of pseudo cycles [7, 9, 18], which shadow long cycles with shorter ones.

For comparison, we also computed the expectation value ⟨x⟩ using the induced map method, the direct cycle expansion, and the stability ordering. In intermittent systems, the presence of a vast number of long-period orbits with extremely small stability eigenvalues leads to extremely slow convergence for both the direct cycle expansion and the stability ordering, regardless of the chosen truncation length ${t}_{{\rm{\max }}}$, and they eventually converge to incorrect values. For the induced map method, once ${t}_{{\rm{\max }}}$ exceeds 10, its accuracy is primarily limited by the approximation degree n_Ind. Since the induced map method and the layered scheme involve different parameters, a direct comparison is not feasible. Instead, we compare them under approximately equivalent computational costs. When ${t}_{{\rm{\max }}}=10$ and n_Ind = 8, the induced map method converges to ⟨x⟩ ≈ 0.363 with an accuracy of about 10^−2.5, slightly worse than the layered scheme (with ${t}_{{\rm{\max }}}=16$, n = 8), which achieves an accuracy of 10^−4.5. Increasing n_Ind to 17 while keeping ${t}_{{\rm{\max }}}=10$ improves the result to ⟨x⟩ ≈ 0.3624 with accuracy 10^−4.5, still slightly below the 10⁻⁵ accuracy attained by the layered scheme with ${t}_{{\rm{\max }}}=16$ and n = 10. Therefore, in intermittent systems the layered scheme significantly outperforms the other three approaches in both accuracy and computational efficiency.

Figures 4(a) and (b) benchmark our layered scheme against the conventional cycle expansion and the stability ordering scheme for a system with four singularities located at the endpoints of the phase interval and the period-2 orbit points as displayed above. To minimize unknown coefficients and enhance stability, we derived analytic relations between singular layer coefficients from equation (5): where ${x}_{{s}_{i}}$ represents the sequence of singularities along the trajectory evolving from the first one at x = 1 (for i = 1), to the later ones at: ${x}_{{s}_{i}}=1,0,{x}_{a},{x}_{b},{x}_{a},{x}_{b}$, each having a specific ${c}_{{s}_{i}}$. Our approach achieves exponential accuracy with increasing t_max and reaches a accuracy level of 10⁻⁶ for ${t}_{\max }=12$, while traditional method and stability ordering shows algebraic convergence [figure 4(b)]. Moreover, it is evident that even with a relatively short truncation length, our scheme can attain a considerably higher numerical accuracy than existing methods. This indicates that, compared to the traditional cycle expansions, which only utilizes information from UPOs, the additional dynamical information introduced here plays an essential role.

For a system with infinitely many singularities [see figure 5(a)], generated by iterating an initial singularity at ${x}_{r}=f(\frac{1}{2})$ through ${f}^{i}({x}_{r})\,(i\in {\mathbb{N}})$, all singular neighborhoods share the same singular exponent $-\frac{1}{2}$. These aperiodic points densely fill the accessible phase space ${ \mathcal M }=[f({x}_{r}),{x}_{r}]$ [see figure 5(b)]. While the measure decomposes into infinitely many singular layers plus a hyperbolic component, solving an infinite-rank equation is computationally intractable. We therefore exploit the self-similar structure of singularities to reduce degrees of freedom and implement recursive relations among the coefficients based on equation (5): where ${f}^{{\prime} }(x)=R-2Rx$ is the stability factor at each singularity of f(x). Besides, a finite truncation of the singularities seems reasonable since the asymptotic decay of singularity peaks upon iteration is clearly observed in figure 5(b), that is, higher-order singularities exhibit exponentially diminishing influence. We implement a simple τ-controlled truncation (here, τ denotes the number of singularities retained sequentially in the computation) and compare the average value ⟨x⟩ computed via a direct cycle expansion (⟨x⟩_Dir), the stability ordering scheme (⟨x⟩_Stab) and our layered approximation (⟨x⟩_Lyr,τ) for τ = 10, 20, 30. figures 6(a) and (b) demonstrate our layered scheme's exponential efficiency. The results of the layered scheme ⟨x⟩_Lyr,τ can attain a considerably high accuracy at a relatively short t_max. By combining the hyperbolic layer's cycle expansion—which supplies the baseline—with the singular layers that provide the necessary corrections, the layered approximation scheme achieves significantly faster convergence overall than other schemes. Furthermore, the truncation of singular layers has an important impact on the convergence. Figure 6(a) reveals that, for small truncation (t_max ≤ 5), the results obtained with τ = 30 are more accurate than those obtained with τ = 10 and τ = 20. Conversely, figure 6(b) shows that, as t_max grows, the convergence of τ = 10 seems to be more steady than that of τ = 30 and τ = 20. A key practical consideration is the trade-off between singularity resolution (large τ) and numerical stability. These trends signal that the retained number τ of singular layers should be tuned to the truncation t_max: at small t_max a larger τ enhances accuracy, whereas for larger t_max an excessive count of singular layers destabilizes the computation and degrades precision, possibly because of the kinks generated by including singular functions in the layered expansion. Practitioners can simply adopt τ = 20 as a robust default without noticeable loss of precision for different ${t}_{{\rm{\max }}}$, adjusting only when high precision at low truncation is required.

The layered approximation scheme, combining additional dynamical information with cycle expansions, ensures rapid convergence and high computation accuracy, resolving the traditional accuracy and efficiency issues to some extent. Decomposing the natural measure into singular and hyperbolic layers clarifies how singularities influence the system's statistical properties, which may play important roles in the analysis of complex non-equilibrium systems. From the perspective of natural measure, the dynamical invariance condition provides extra corrections beyond the UPOs. This opens a promising path towards wider application of cycle expansions and high-precision computation for chaotic systems.

In this work, we employed the layered scheme to study dynamical systems with distinct singularity characteristics: (i) Intermittent systems, where the natural measure possesses a singularity at the origin; (ii) multi-singularity systems, where singularities arise both on periodic orbits and at phase-space boundaries; (iii) systems exhibiting infinitely many singularities in the natural measure. In these numerical examples, the numerical accuracy and convergence rate of the layered scheme are superior to those of existing schemes and hence relieve a few well-known limitations of traditional cycle expansions and existing schemes.

Nevertheless, the layered scheme still has its own limitations. How to achieve an elegant trade-off between numerical accuracy and stability within this method still calls for deeper investigation; at present we can only proceed by trial-and-error or empirical tuning. While the current work focuses on one-dimensional maps for clarity and benchmarking, its application to more complex systems still requires further investigation. The conceptual framework of decomposing the invariant measure into singular and hyperbolic components remains conceptually applicable. In higher dimensions, singularities often manifest on unstable manifolds or fractal basin boundaries [33, 34]. The challenge lies in analytically characterizing these complex singularity structures and constructing a suitable hyperbolic conjugate system. This is a promising direction for future work, potentially building on frameworks for cycle expansions in higher-dimensional systems [5, 35–37]. However, the convergence issues of cycle expansion have not yet been fully resolved, and balancing mathematical rigor with physical operability remains an open question. The Layered scheme may provide useful insights for the readers and offer a potential approach to tackling more realistic challenges in the future.

This work was supported by the National Natural Science Foundation of China (No. 12375030), by BUPT Excellent PhD Students Foundation (CX2023210), and also by the Key Program of National Natural Science Foundation of China (No. 92067202).