Rigorous and effective numerics for Liverani-Saussol-Vaienti maps

Alexey Korepanov; YuTong Wei; Caroline Wormell

Source-linked AI summary

Rigorous and effective numerics for Liverani-Saussol-Vaienti maps

Alexey Korepanov, YuTong Wei, Caroline Wormell

arXiv:2610.01879v1Oct 1, 2026math.DScs.NAmath.NA

TL;DR

Accurate estimation of LSV-map invariant densities is challenging because nonuniform expansion creates an infinite, slowly decaying induced-operator sum. The paper combines an Abel function with Euler–Maclaurin summation and finite-dimensional operator approximation. Its method has polynomial accuracy-cost growth and establishes analyticity of the invariant density for all γ>0.

Problem
Rigorous, accurate estimation of LSV-map invariant densities is difficult because nonuniform expansion makes the induced transfer operator an infinite, slowly decaying sum over iterates.
Source: §1 · p. 1
Method
The approach uses an asymptotically computed Abel function, Euler–Maclaurin summation, and a finite Chebyshev matrix approximation of the transfer operator.
Source: §1 · p. 2
Results
The method’s computational effort grows polynomially with digits of accuracy, and the invariant density is real-analytic on (0, 1] for every γ>0.
Source: §1 · p. 2
Takeaways & Limitations
Numerical monotonicity checks for h′/h and h′′/h support a relatively simple proof of optimal memory-loss rates for finite and infinite invariant measures.
Source: §1 · p. 2
The average-return-time integral is treated under the assumption γ<1 so that it is finite.
Source: §5.1 · p. 20

Abstract

from arXiv · show

We make highly accurate numerical estimates of basic properties of the invariant measure for Liverani-Saussol-Vaienti intermittent maps, an archetypal slowly mixing dynamical system. We solve the challenge of precise and efficient estimation for a map with non-uniform expansion, where the transfer operator does not have a spectral gap. We do this using an Abel function, which solves the neutral dynamics, and which we can compute accurately via an asymptotic expansion. Our work covers both finite and sigma-finite invariant measure cases. In particular, we obtain the first practical estimates for the parameter far from zero, including the sigma-finite case. This opens the door to much deeper numerical study of intermittent dynamics than was previously possible. A driving motivation for this work is the upcoming proof of optimal rates of memory loss which works equally for finite and infinite invariant measure cases. No prior results on rates of memory loss for dynamical systems with infinite invariant measure exist, with a single exception of a recent paper by I.Chevyrev and A.K. on null recurrent Markov chains.

1 Introduction

The paper develops rigorous, accurate numerical estimates of invariant densities and related properties for LSV maps despite nonuniform expansion and slowly computed induced-operator sums. Its Abel-function-based approach supports estimates across a wide range of parameters and motivates numerical verification of conditions for memory-loss results.

1 Introduction

The goal is rigorous, accurate estimation of h, a challenge because nonuniform expansion makes the induced transfer operator an infinite, slowly decaying sum over iterates.For each γ>0, the invariant measure is finite when γ<1 and infinite when γ≥1; no closed forms are known for it or most related statistical quantities.Source: §1 · p. 1
An Abel function resolves the neutral dynamics; Euler–Maclaurin summation and a finite Chebyshev matrix then approximate the transfer operator and invariant density.The Abel function has a quickly converging asymptotic expansion for computing orbits near zero.Source: §1 · p. 2
The method’s computational effort grows polynomially with digits of accuracy, unlike earlier exponential-cost methods that practically limited estimates to γ ≲1.The approach also requires effort that does not vary massively with γ.Source: §1 · p. 2
Numerical checks of h′/h and h′′/h monotonicity on (0, 1] support a relatively simple proof of optimal memory-loss rates for finite and infinite invariant measures.The memory-loss results themselves are reserved for a follow-up paper in preparation.Source: §1 · p. 2
The invariant density of the LSV map is real-analytic on (0, 1] for every γ>0.The paper also notes that its results generalize to maps f(x)=x f̂(x^γ) near zero when f̂ is analytic.Source: §1 · pp. 2–3

Code and examples

Examples report a computation taking a couple of days and a low-precision computation taking seconds on a laptop or minutes in a web browser; code and a web demo are public.The examples include average return time and Lyapunov-exponent computations.Source: Code and examples · p. 3

2 Abel function

The section constructs an Abel function for the left branch and proves it is analytic near zero and can be computed accurately and efficiently. It establishes the functional relation, convergence of the construction, and bounds needed for later work.

2 Abel function

The Abel function for the left branch is analytic near zero and accurately and efficiently computable, with approximants chosen to satisfy its functional relation to increasing order.The function is constructed using the coordinate z = x^γ and an asymptotic expansion.Source: §2 · pp. 3–4
The defining limit for the Abel function is analytic on the petal and converges uniformly on compact subsets, independently of the approximant order.The proof bounds the residual terms by a summable sequence along inverse iterates.Source: §2 · p. 5

2.1 Estimation of M

A positive-coefficient majorant bounds the functional-equation defect on the petal, supporting the explicit estimation of the constant M.The coefficients of the majorant dominate those of the defect, and its bound increases with |z|.Source: §2.1 · p. 6

2.2 t-domain and derivatives

In the t-coordinate, the analysis bounds derivatives of the Abel function and establishes that its derivative is negative on the positive real axis.The Abel function also extends to points that enter the petal under inverse iteration.Source: §2.2 · pp. 6–7

3 Sums for induced transfer operator

This section develops estimates for sums arising in the transfer operator of the first return map to [1/2, 1]. It applies an Euler–Maclaurin formula and derives bounds for bounded and specific observables.

3 Sums for induced transfer operator

The sums considered arise in the transfer operator of the first return map to [1/2, 1], and their estimates reduce computation for large x to small x.The section applies its estimates to particular choices of the observable φ(x, k).Source: §3 · p. 8, §3.4 · p. 12

3.1 Euler-Maclaurin formula

The adapted Euler–Maclaurin result bounds the sum approximation remainder under analyticity and uniform disk-bound assumptions.The paper takes the abstract result with a minor modification from [VW25, Section 5].Source: §3.1 · p. 8

3.2 Preliminary bounds

Preliminary bounds show that the relevant shifted points remain in the image of the petal under the Abel function, enabling control of the integrand there.The domain argument follows paths in the image and its inverse under the Abel function.Source: §3.2 · p. 9

3.3 Bounded observables

For bounded observables, the section bounds the analytic function required by the Euler–Maclaurin estimate and combines that bound with the remainder estimate.The setup includes observables such as Chebyshev polynomials evaluated along inverse iterates.Source: §3.3 · pp. 10–11

3.4.1 Chebyshev polynomials

For Chebyshev-polynomial observables, the section uses a maximum bound and recursively computable integrals in its numerical estimates.The maximum of φ_m on [-δ, δ] is attained at -δ.Source: §3.4 · p. 12, §3.4.1 · p. 12

3.4.2 First return time

For first-return-time observables, approximating A with its asymptotic approximant and using |φ_m| ≤ 1 yields a bound on the integral error.The observable is φ(x, k) = (k + K)φ_m(x^k).Source: §3.4.2 · pp. 12–13

3.4.3 Lyapunov exponent

For the Lyapunov-exponent observable, the section expands the logarithmic term so the resulting integrals of powers times φ_m have closed forms.The observable is φ_m(x^k) log J_k(x).Source: §3.4.3 · pp. 13–14

4 Invariant density on [1/2, 1]

The section establishes a controlled Chebyshev representation of the induced invariant density on [1/2, 1], using complex-domain bounds and a Hardy-space operator framework.

4 Invariant density on [1/2, 1]

The induced invariant density is analytic on ellipse A and admits a finite Chebyshev approximation with an explicit error bound.Source: §4 · p. 14

4.1 Bernstein ellipses

Every inverse branch of the induced map sends ellipse C into A, supporting the analytic construction of the density approximation.Source: §4.1 · pp. 14–15

4.2 Bounds on sums in A

The sums considered in A are finite when γ < 1.Source: §4.2 · p. 16

4.3 Bounds in Hardy spaces

On Hardy space HA, L has a maximal eigenvalue 1 of multiplicity one and a spectral gap, with a unique normalized invariant density.Source: §4.3 · p. 17
The approximation ȟ has an explicit error bound controlled by the resolvent norm, the operator approximation error, and ∥ȟ∥A.Source: §4.3 · p. 20

5 Return time and Lyapunov exponent

The section develops rigorous numerical estimates for the mean return time and Lyapunov exponent using the induced density and transfer-operator approximations.

5 Return time and Lyapunov exponent

The section's principal applications are estimating the mean return time to [1/2, 1] and the Lyapunov exponent.Source: §5 · p. 20

5.1 Return time

The mean return time is computed from the first return time τ and invariant density h, with γ < 1 required for the defining integral to be finite.Source: §5.1 · p. 20

5.2 Lyapunov exponent

Estimating Λ is identified as the primary difficulty in computing the Lyapunov exponent, and the section assembles bounds for its numerical approximation.Source: §5.2 · p. 21

6 Invariant density on (0, 1]

The section extends the density estimate from [1/2, 1] to (0, 1] by combining Euler–Maclaurin evaluation with small- and moderate-x approximations and explicit error control.

6 Invariant density on (0, 1]

The Euler–Maclaurin formula and the approximation on [1/2, 1] are used to compute h(x) accurately for small x.Source: §6 · p. 22
For not-so-small x, the section instead develops a less accurate but rigorous approximation using a crude truncated Abel function.The truncation uses a relatively small iterate count, such as 8, to allow a larger petal.Source: §6 · p. 22
For x ≤ ẍ*, the density derivatives are approximated by terms involving M0 and M1, with explicit bounds on the error terms.Source: §6 · pp. 23–24
The stated derivative bounds require ẍ* to be sufficiently small and M0 > δ0.Source: §6 · p. 24

7 Conditions for memory loss

The section sets up a rigorous-numerics verification of derivative conditions used for memory-loss results, extending bounds from near zero across the intermediate interval.

7 Conditions for memory loss

The derivative properties needed for memory-loss results are not known analytically, so the authors verify them using rigorous numerics.Source: §7 · p. 25
Recursive evaluation of h on the intermediate interval (ẍ*, 1/2) loses too much precision to be practical.Source: §7 · p. 25
The proposed verification checks derivative signs on [1/2, 1] and propagates bounds from a known interval to its image under f.Source: §7 · p. 25
The propagation lemma assumes positive lower bounds for derivatives through order three and an upper bound on −h′ on an interval within (0, 1/2].It then supplies lower bounds for those derivatives on the image interval.Source: §7 · p. 25

A Basic lemmas

The basic lemmas establish coefficient recurrences and a condition involving the complex variable t, with proof steps reducing an inequality to H(R) > 0. A further lemma considers power-series coefficients under a summability condition.

A Basic lemmas

For t with real part at least b(γ + 1), Lemma A.1 states a result for g(t) = t(1 + b/t)−γ, where b = 2γ.The supplied passage gives the lemma's setup but not its full conclusion.Source: A Basic lemmas · p. 26
The proof for t = u + iv reduces a condition to H(R) > 0 and establishes H′(R) ≥ 0 under the stated bounds.The argument uses polar coordinates for 1 + b/t and bounds involving R sin φ and φ.Source: A Basic lemmas · p. 26
A further lemma focuses on the case where a_k = 0 for all sufficiently large k and assumes a weighted coefficient sum is less than a_0.The condition is stated as Σ|a_k|δ^k < a_0 for some δ > 0.Source: A Basic lemmas · p. 26
Equating formal power-series coefficients yields the recurrence for c_k, and the related coefficients satisfy c̃_k ≥ |c_k|.The proof assumes a_0 = 1.Source: A Basic lemmas · p. 27

B Computation of Integrals

The appendix describes algorithms for computing definite integrals of scaled Chebyshev polynomials over [0, x].

B Computation of Integrals

The algorithms compute integrals over [0, x] of Chebyshev polynomials scaled to [0, 1], φ_m(u) = T_m(2u − 1).The appendix describes these algorithms as those used in the authors' code.Source: B Computation of Integrals · p. 27

B.1 Fractional power integrals

The fractional-power integral computation uses a recurrence derived from a Chebyshev-polynomial relation and integration by parts.

B.1 Fractional power integrals

The computation concerns integrals of the form introduced for this subsection, with the recurrence enabling evaluation for m ≥ 2.The displayed integral expression and the explicitly computed cases are not fully specified in the supplied passages.Source: §B.1 · p. 27
For β > −1, integrating the relation between φ_m and its neighboring polynomials by parts yields a recurrence for the fractional-power integrals.The recurrence supports computation for m ≥ 2 after explicitly computing initial cases.Source: §B.1 · p. 27

B.2 Logarithmic integrals

The logarithmic-integral computation differentiates the fractional-power recurrence with respect to β and evaluates it at β = 0.

B.2 Logarithmic integrals

Differentiating the recurrence for I_m,β and evaluating at β = 0 gives a recurrence coupling L_m with the non-logarithmic integrals I_m,0.The subsection applies a technique similar to the fractional-power integral computation.Source: §B.2 · p. 27
The initial logarithmic integrals are given explicitly as L_0 = x log x − x and L_1 = (x^2 − x) log x − 1.Source: §B.2 · p. 27