Struct that encapsulates the additional quantities needed on the coarse basis for a two-grid estimate, or on the (only) basis for a one-grid estimate. It is meant as an intermediate quantity that can be saved on the disk to avoid recomputing Q all the times.

Struct that encapsulates all the quantities computed from the fine basis that are needed in the two-grid estimate. It is meant as an intermediate quantity that can be saved on the disk to avoid recomputing Q all the times

Compute FineGridQuantities and CoarseGridQuantities, given a function f(n) that computes B, D, Q = f(n)

Compute FineGridQuantities, given a function f(n) that computes B, D, Q = f(n)

Bounds rigorously the distance of w from the fixed point of Q (normalized with integral = 1), using a vector of bounds norms[k] ≥ ||Qh^k|{U_h^0}||. If ε₁ and normQ are given, then Q can be omitted

distance_from_invariant_residual(B::Basis, D::Dynamic, Q, w, norms;
    ε₁ = residualbound(B, weak_norm(B), Q, w),
    ε₂ = mag(integral_covector(B) * w - 1),
    dfly_coefficients = dfly(strong_norm(B), aux_norm(B), D),
    defect_coefficients = projection_defect_coefficients(B, D; dfly_coefficients))

A posteriori (residual-based) rigorous upper bound for the weak-norm distance $\|h - w\|_w$ between the invariant density $h$ of the abstract transfer operator $L$ of the dynamic D and a computed candidate w in the approximating space of B.

This is the "exchanged" version of [Galatolo–Monge–Nisoli–Poloni, Chaos Solitons & Fractals 170 (2023) 113329, Theorem 3.4] anticipated in Remark 3.9 of that paper: there, the error $v = h - w$ is expanded in powers of the discretized operator $Q_h$, and one pays the a priori strong-norm bound $\|h\|_s \le B/(1-A)$ (Corollary 2.6) for the unknown density — that is distance_from_invariant. Here the error is instead expanded in powers of the abstract operator $L$, driven by the computed residual of the candidate; Remark 3.9 notes this requires summing $\|L^k r\|_w$, "a difficult task" a priori — the closed-form two-norm bound below performs that summation using only the computed norms $C_k$ of the discretized operator, the one-step Lasota–Yorke inequality, and the discretization-defect constants. It is sharper than the a priori bound whenever $\|w\|_s \ll B/(1-A)$, i.e. precisely when the Lasota–Yorke constants are poor ($A$ near $1$, small branches) — the candidate's strong norm is computed, not estimated.

Framework and assumptions (notation of the cited paper)

• $L$ preserves the integral, $\|L\|_{L^1}\le 1$, and satisfies the one-step Lasota–Yorke inequality $\|Lf\|_s \le A\|f\|_s + B\|f\|_{L^1}$ with $A < 1$ (dfly_coefficients, cf. dfly); the norms satisfy $\|\cdot\|_{L^1} \le \|\cdot\|_w$ (Assumption 2.1(4) of the paper).
• Q is the rigorously assembled, integral-preserving discretized operator $Q_h$ of a compatible discretization; norms[k] $\ge \|Q_h^k|_{V_h^0}\|_w$ are certified (cf. powernormbounds, finepowernormbounds); $K_h$ = weak_projection_error.
• defect_coefficients $= (c_s, c_w)$ satisfy $\|(L-Q_h)v\|_w \le K_h(c_s\|v\|_s + c_w\|v\|_{L^1})$ (see projection_defect_coefficients).
• $L$ has a unique invariant probability density $h$; uniqueness identifies the Neumann series below with $\bar w - h$ (any integral-zero fixed point of $L$ vanishes by ergodic decomposition). In practice this is certified from the same norms by the small-matrix method (convergencerateabstract, Galatolo–Nisoli–Saussol).
• The strong and $L^1$ norms of the candidate are bounded rigorously via normbound(B, strong_norm(B), w) and normbound(B, aux_norm(B), w); the basis must provide these methods.

Derivation

Let $\bar w = w/\int w$ (normalization defect $\varepsilon_2$ accounted at the end) and $r = L\bar w - \bar w$ the abstract residual, $i(r) = 0$. Then

\[\bar w - h = -(I-L)^{-1}\big|_{V^0}\, r = -\sum_{k\ge0} L^k r, \qquad \|h - \bar w\|_w \le \sum_{k\ge0} W_k,\]

with $W_k := \|L^k r\|_w$, $S_k := \|L^k r\|_s$.

(1) Both norms of the residual are computable. Writing $r = (L - Q_h)\bar w + (Q_h\bar w - \bar w)$, the second term is the computed eigen-residual ($\le \varepsilon_1/(1-\varepsilon_2)$) and the first pays the defect of the explicit candidate:

\[W_0 \le \frac{\varepsilon_1 + K_h (c_s \|w\|_s + c_w \|w\|_{L^1})}{1-\varepsilon_2}, \qquad S_0 \le \frac{(1+A)\|w\|_s + B\|w\|_{L^1}}{1-\varepsilon_2},\]

the latter from one Lasota–Yorke step applied to w and the triangle inequality.

(2) Two-channel recursion. Telescoping $L^k = Q_h^k + \sum_{j<k} Q_h^j (L - Q_h) L^{k-1-j}$ (all vectors have zero integral, so the restricted norms apply) and iterating the Lasota–Yorke inequality:

\[S_k \le A^k S_0 + B\sum_{i<k} A^i W_{k-1-i},\qquad W_k \le C_k W_0 + K_h \sum_{j<k} C_j\,\big(c_s S_{k-1-j} + c_w W_{k-1-j}\big).\]

(3) Sound closed-form summation. Summing over $k \ge 0$ (Fubini for nonnegative series), with $\Sigma_C := \sum_{k\ge0} C_k$ bounded by infinite_sum_norms:

\[\Sigma_S \le \frac{S_0}{1-A} + \frac{B}{1-A}\Sigma_W,\qquad \Sigma_W \le \Sigma_C W_0 + K_h\Sigma_C\big(c_s\Sigma_S + c_w\Sigma_W\big),\]

whence, provided the denominator below is positive (checked; an error is thrown otherwise),

\[\|h - \bar w\|_w \;\le\; \Sigma_W \;\le\; \frac{\Sigma_C W_0 + K_h \Sigma_C\, c_s\, S_0/(1-A)} {1 - K_h \Sigma_C \big(c_s\, B/(1-A) + c_w\big)} .\]

The same bound applies to every partial sum, so the Neumann series converges absolutely and the identity above is justified. Finally $\|h - w\|_w \le \Sigma_W + \tfrac{\varepsilon_2}{1-\varepsilon_2}\|w\|_w$.

All arithmetic is carried out in interval arithmetic; the returned value is a rigorous Float64 upper bound.

Uses power norm bounds already computed for a coarse operator to estimate the same norms for a finer operator

Return a numerical approximation to the (hopefully unique) invariant vector of the dynamic with discretized operator Q.

The vector is normalized so that integral_covector(B)*w ≈ 1

Compute a one-grid error estimate.

The first return argument is the error, the second is the time breakdown according to ["dfly", "assembling", "eigen", "norms", "error"]. (The sum of that vector is the total time taken)

Uses different strategies to compute power norm bounds.

If specified, m norms of powers are estimated computationally, and then m_extend norms are obtained with a cheaper refinement process. Otherwise these numbers are selected automatically.

A vector of length mextend is returned, such that norms[k] ≥ ||Qh^k|{Uh^0}||

Computes bounds for norms of powers, taking (optionally) minimum values for the number of norms to compute

projection_defect_coefficients(B::Basis, D::Dynamic;
    dfly_coefficients = dfly(strong_norm(B), aux_norm(B), D),
    normL = bound_weak_norm_abstract(B, D; dfly_coefficients = dfly_coefficients))

Certified constants $(c_s, c_w)$ such that the discretization defect of the (integral-preserving) discretized operator $Q_h$ satisfies, for every $v$ in the strong space,

\[\|(L - Q_h)\,v\|_w \;\le\; K_h \,\big( c_s \|v\|_s + c_w \|v\|_{L^1} \big), \qquad K_h = \texttt{weak\_projection\_error}(B).\]

Generic method, valid for every compatible discretization in the sense of [Galatolo–Monge–Nisoli–Poloni, Chaos Solitons & Fractals 170 (2023) 113329, Definition 2.7]: by Lemma 3.5 of that paper, $\|(Q_h - L)f\|_w \le 2K_h(\|L\|_w \|f\|_s + \|Lf\|_s)$, and the one-step Lasota–Yorke inequality (5) gives $\|Lf\|_s \le A\|f\|_s + B\|f\|_{L^1}$, whence

\[(c_s, c_w) = \big(2(\|L\|_w + A),\; 2B\big).\]

Sharper basis-specific methods may be provided (see the Ulam method, where $(c_s, c_w) = (1+A, B)$).

Return an upper bound to Q_h*w - w in the given norm

Compute a two-grid error estimate.

The first return argument is the error, the second is the time breakdown according to ["dfly", "coarse", "assembling", "eigen", "norms", "error"]. (The sum of that vector is the total time taken)

analytic_dfly(strong::Aη, η′, C, K) -> (A, B)

Lasota–Yorke constants for the Fourier analytic norm $\|f\|_{A_η} = \sum_k |\hat c_k| e^{2πη|k|}$, with auxiliary norm $L^1(dx)$:

\[\|Lf\|_{A_η} \le A\,\|f\|_{A_η} + B\,\|f\|_{L^1}.\]

η′ > η is the enlarged strip certified by annulus_expansion and C bounds $\|L\|_{A_η \to A_{η′}}$. Splitting the Fourier series at K:

\[A(K) = C\,\frac{2q^{K+1}}{1-q},\quad q = e^{-2π(η′-η)}, \qquad B(K) = \sum_{|k| \le K} e^{2πη|k|} .\]

The B term uses $|\hat c_k(Lf)| \le \|Lf\|_{L^1} \le \|f\|_{L^1}$: the coefficient bound is uniform in k, and the transfer operator is an $L^1$ contraction, so no constant is lost.

analytic_dfly(strong::Eρ, ρ′, C, K; L1μ_bound = 1.0) -> (A, B)

Lasota–Yorke constants for the Chebyshev analytic norm $\|f\|_{E_ρ} = \sum_k |\hat b_k| ρ^k$, with auxiliary norm $L^1(μ)$, μ the arcsine measure:

\[\|Lf\|_{E_ρ} \le A\,\|f\|_{E_ρ} + B\,\|f\|_{L^1(μ)} , \qquad A(K) = C\,\frac{2q^{K+1}}{1-q},\; q = ρ/ρ', \qquad B(K) = 2 D_μ \sum_{k \le K} ρ^k .\]

ρ′ > ρ is the enlarged ellipse certified by bernstein_expansion and C bounds $\|L\|_{E_ρ \to E_{ρ'}}$.

analytic_dfly_choose_K(strong, gain, C; target_A = 0.5, max_K = 4096, kwargs...)

Smallest K for which analytic_dfly achieves A ≤ target_A, returned as (K, A, B). Since A(K) decreases geometrically and B(K) grows geometrically, this is the cheapest usable inequality: any smaller K fails the target, any larger one inflates B for nothing.

Returns nothing if the target is not reachable within max_K.

analytic_dfly_degenerate(strong, gain, C₂) -> (A, 0.0)

The degenerate Lasota–Yorke $\|Lf\|_{A_η} \le A\,\|f\|_{A_η}$, with no auxiliary term — that is, A is simply the continuity constant of L on the analytic space, valid for every f, not on any subspace.

C₂ is the bound between the two domains, $\|L\|_{A_η \to A_{η'}}$, and the enlargement is converted into a constant by splitting the weight at the gain δ:

\[\|Lf\|_{A_η} = \sum_k e^{-2πkδ}\,e^{2πkη'}|\hat c_k(Lf)| \;\le\; \Big(\sum_k e^{-2πkδ}\Big)\,\|Lf\|_{A_{η'}} \;\le\; G(δ)\,C_2\,\|f\|_{A_η},\]

the geometric sum converging exactly because the neighbourhood is enlarged. For Eρ the weight is ρ^k and the ratio is ρ/ρ'.

This is the single-space setting of Nisoli, Certified spectral approximation of transfer operators and the Gauss map, arXiv:2602.19435, where the analytic (Hardy-space) case needs only $\|L\|_{B \to B} \le C$ together with the truncation bound coming from the domain enlargement — the strong–weak DFLY scale of that paper's Appendix A being a separate setting. Compactness, not a contraction factor, is what drives the spectral certification, so A here need not be less than 1.

analytic_transfer_bound(min_derivative, nbranches) -> C

$C$ with $\|Lf\|_{∞,\,\text{nbhd}} \le C\,\|f\|_{∞,\,\text{nbhd}}$, for the transfer operator $Lf(z) = \sum_k f(g_k(z))\,g_k'(z)$.

Since $|g_k'| = 1/|T'\circ g_k|$ and the inverse branches land in the neighbourhood where min_derivative was certified,

\[C \;\le\; \frac{\#\text{branches}}{\min |T'|} .\]

min_derivative comes from min_modulus_on_ellipse or min_modulus_on_circle — i.e. from the complex neighbourhood, never from the assembled matrix.

analytic_truncation_error(strong, gain, C₂, K) -> ε_K

$ε_K \ge \|(I - Π_K)L\|_{\mathcal B \to \mathcal B}$ on the analytic space $\mathcal B$, i.e. the truncation number of Assumption 2.1(3).

One line, from the domain gain alone. For f in the unit ball of $E_ρ$, $Lf$ lies in $A(E_{ρ'})$ with $\|Lf\|_{E_{ρ'}} \le C_2$; writing $Lf = \sum_j b_j T_j$,

\[\|(I-Π_K)Lf\|_{E_ρ} = \sum_{j>K}|b_j|ρ^j = \sum_{j>K}|b_j|ρ'^j\Big(\frac{ρ}{ρ'}\Big)^{j} \le C_2\,q^{K+1}, \qquad q = ρ/ρ' < 1 .\]

The Fourier case is identical with $q = e^{-2π(η'-η)}$.

annulus_expansion(f, η; n = 1024) -> (η_in, η_out)

Certified enlargement of the strip $|\mathrm{Im}\,x| \le η$ under a circle map f, the Fourier counterpart of bernstein_expansion.

Under $x \mapsto z = e^{2πix}$ the strip becomes the annulus $e^{-2πη} \le |z| \le e^{2πη}$, so the image of each boundary circle is enclosed and its extreme modulus recorded:

\[η_{out} = \frac{\log \min_{|z| = e^{2πη}} |f(z)|}{2π}, \qquad η_{in} = \frac{-\log \max_{|z| = e^{-2πη}} |f(z)|}{2π}.\]

f acts on the z variable and must accept a Complex{Interval}. The strip is enlarged when both exceed η; min(η_in, η_out) is then the certified η′ for analytic_dfly.

bernstein_geometry(strong::Eρ, D::PwMap; n = 1024) -> (ρ′, C₂)

The two numbers the analytic theory rests on, computed once from the map:

• ρ′ > ρ — the enlarged Bernstein parameter, certified by bernstein_expansion applied to each branch in the GLOBAL chart $t = 2x-1$. It means the inverse branches contract: $g_k(E_{ρ'}) ⊆ E_ρ$.
• C₂ = Σ_k 1/\min_{E_ρ}|T_k'| — a bound for $\|L\|_{A(E_ρ) → A(E_{ρ'})}$.

These are exactly the quantities dfly(::Eρ, ::Type{L1}, ::PwMap) computes internally and then discards; exposing them avoids recomputing.

dfly(norm::Aη, ::Type{L1}, D::PwMap; C₂, n = 1024) -> (A, 0.0)
dfly(norm::Eρ, ::Type{L1}, D::PwMap; C₂, n = 1024) -> (A, 0.0)

Degenerate Lasota–Yorke for the analytic norms: A is the continuity constant of L on the analytic space, and the auxiliary constant is 0.

Same call shape as every other dfly, so powernormbounds, refine_norms_of_powers and the coarse–fine workflow run unchanged once the basis fixes its strong and weak norms. Because B = 0 the auxiliary norm is multiplied by zero and plays no part; the weak norm is free to be whatever the approximation error is measured in, which for these bases is L2.

The geometry needs only forward enclosures — the branches of a PwMap evaluate at Complex{Interval} directly, so the enlarged neighbourhood is certified by strip_expansion / bernstein_expansion applied to each branch reparametrized onto its own domain. No complex derivative is involved.

C₂ bounds L between the two domains, $\|L\|_{A_η \to A_{η'}}$; the Hölder factor G(δ) from the enlargement then gives A = G(δ)C₂ (see analytic_dfly_degenerate).

C₂ is the one ingredient the forward enclosure does not determine on its own: the transfer operator carries the weight $|g_k'| = 1/|T_k'\circ g_k|$, so by default it is computed as $\sum_k 1/\min|T_k'|$ with the minimum taken over the enclosed neighbourhood — valid on the closed neighbourhood by the minimum modulus principle, $T_k'$ being zero-free for an expanding map. Pass C₂ explicitly to override it with a sharper estimate for your operator.

(A, B) = dfly(strongnorm, auxnorm, dynamic)

Constants (A, B) such that ||Lf||s ≦ A||f||s + B||f||_aux

invariant_measure_strong_norm_bound(B, D; dfly_coefficients)

Strong-norm bound on the invariant density.

The classical DFLY route is $B/(1-A)$, which needs A < 1 and B > 0. In the degenerate analytic case B = 0 and A is a continuity constant, and the bound is simply A — the density lies in the image of L, so the continuity constant already controls it.

min_modulus_on_circle(g, η; n = 1024)

Lower bound on |g| over both circles |z| = e^{±2πη}, the strip counterpart of min_modulus_on_ellipse.

min_modulus_on_ellipse(g, ρ; n = 1024)

Rigorous lower bound on $\min_{∂E_ρ} |g|$, by the same boundary covering as bernstein_expansion.

If g has no zero in the closed ellipse — which is the case for T' of an expanding map — the minimum modulus principle makes this a bound on the whole of $E_ρ$, not just its boundary.

single_space_resolvent_lift(R_K, ε_K) -> R

Corollary 2.3 of the Gauss-map paper: if $α := ε_K R(z,L_K) < 1$ then $z ∈ ρ(L)$ and

\[R(z, L) \;\le\; \frac{R(z, L_K)}{1 - ε_K R(z, L_K)} .\]

Returns Inf when α ≥ 1.

Both arguments must be measured in the SAME norm — for the analytic bases, the $ℓ^2$ Bernstein norm of bernstein_l2_resolvent_bound, which is a diagonal reweighting and so is directly accessible to the verified SVD.

The result is a certificate on the ABSTRACT operator L. It is therefore a property of the map, not of the discretization: computed once at a small K it may be reused at any basis size and any working precision, and never needs recomputing.

strip_expansion(f, η; n = 1024) -> η′

min of the two values from annulus_expansion: the largest strip half-width certified to be reached from η.

strip_geometry(strong::Aη, D::PwMap; n = 1024) -> (η′, C₂)

Fourier counterpart of bernstein_geometry: the enlarged strip half-width and $\|L\|_{A_η → A_{η'}}$.