AverageZero{B<:Basis}

Yield a basis of the space of average zero vectors

Replacement of DualComposedWithDynamic.

aux_normalized_projection_error(B::Basis)

Return a constant Eh (typically scales as h ~ 1/n) such that

$|||P_h f|||\leq |||f|||+ Eh * ||f||_s$

Must be rounded up correctly!

aux_weak_bound(B::Basis)

Return a constant $M₂$ such that for a vector $v ∈ Uₕ$

$|||v|||\leq M_2||v||$

Must be rounded up correctly!

bound_linalg_norm_L1_from_weak(B::Basis)

Return a constant $A$ such that for a vector $v ∈ Uₕ$

$||v||_1\leq A||v||$

Must be rounded up correctly!

bound_linalg_norm_L∞_from_weak(B::Basis)

Return a constant $A$ such that for a vector $v ∈ Uₕ$ ||v||_\infty \leq A||v|| Must be rounded up correctly!

bound_weak_norm_abstract(B::Basis, D=nothing; dfly_coefficients=nothing)

Returns an a priori bound on the weak norm of the abstract operator $L$

bound_weak_norm_from_linalg_norm(B::Basis)

Return constants W₁, W₂ such that for a vector $v ∈ Uₕ$

$||v||\leq W_1||v||_1+W_2||v||_{\infty}$

Must be rounded up correctly!

evaluate(B::Basis, i, x)

Evaluate the i-th basis element at x

evaluate_integral(B::Basis, i; T = Float64)

Value of the integral on [0,1] of the i-th basis element

integral(B::Basis, v; T = Float64)

Return the integral of the function with coefficients v in the basis B

integral_covector(B::Basis)

Return a covector that represents the integral in the basis B

invariant_measure_strong_norm_bound(B::Basis, D::Dynamic)

Bounds $||u||_s$, where $u$ is the invariant measure normalized with $i(u)=1$.

is_integral_preserving(B::Basis)

Integral-preserving discretizations may specialize this to "true"

is_refinement(Bfine::Basis, Bcoarse::Basis)

Check if Bfine is a refinement of Bcoarse

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> Bfine = Ulam(2048)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 2049))

julia> is_refinement(Bfine, B)
true

julia> Bfine = Ulam(2049)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 2050))

julia> is_refinement(Bfine, B)
false

one_vector(B::Basis)

Vector that represents the function 1 in the basis B

projection(B::Basis, f::Function; kwargs...)

Project a function f onto the basis B. The result type and available keyword arguments are basis-dependent. Currently implemented via the TaylorModels extension for the Ulam basis.

strong_norm(B::Basis)

Return the type of the strong norm of the basis

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> strong_norm(B)
TotalVariation

strong_weak_bound(B::Basis)

Return a constant $M₁n$ such that for a vector $v ∈ Uₕ$

$||v||_s\leq M1n*||v||$

Must be rounded up correctly!

Return constants $S₁, S₂$ such that for a vector $v ∈ Uₕ$

$||v||\leq S_1||v||_s+S_2|||v|||$

Must be rounded up correctly!

weak_norm(B::Basis)

Return the type of the weak norm of the basis

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> weak_norm(B)
L1

weak_projection_error(B::Basis)

Return a constant Kh (typically scales as h ~ 1/n) such that

$||P_h f-f||\leq Kh ||f||_s$

Must be rounded up correctly! This function is not exported explictly but is used in all the estimates.

Example

julia> using RigorousInvariantMeasures;

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> RigorousInvariantMeasures.weak_projection_error(B)
0.00048828125

Ulam

Ulam basis on [0,1] associated to the partition $p = \{x_0 = 0, x_1, \ldots, x_n=1\}$

Ulam(n::Integer)

Equispaced Ulam basis on [0,1] of size n

p1 * p2  for two `ProjectedFunction{<:Ulam}` on the same basis

Componentwise multiplication of the cell-value vectors. For Ulam each $p_i.v$ is the cell-value vector of the piecewise-constant reconstruction; their pointwise product is again piecewise-constant on the same partition with cell value $p_1.v[i] \cdot p_2.v[i]$.

Bounds combine via Hölder:

\[\|fg - φ_{p_1}φ_{p_2}\|_{L^1} \;\leq\; \|f\|_{L^∞}\,\|g - φ_{p_2}\|_{L^1} + \|φ_{p_2}\|_{L^∞}\,\|f - φ_{p_1}\|_{L^1}\]

so

weak_dual_bound = p1.weak_dual_bound * p2.weak_dual_bound
proj_error =
    p1.weak_dual_bound * p2.proj_error
    + p2.weak_dual_bound * p1.proj_error

Base.getindex(B::Ulam, i::Int)

Returns the i-th element of the Ulam basis as a function.

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(16)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 17))

julia> B[1](1/32)
1

julia> B[2](1/32)
0

iterate(S::AverageZero{Ulam{T}}, state = 1) where{T}

Return the elements of the basis of average $0$ functions in the Ulam Basis

Given a preimage of an interval I_i, this iterator returns its relative intersection with all the elements of the Ulam basis that have nonzero intersection with it

Base.length(B::Ulam)

Returns the size of the Ulam basis (the size of the underlying vector -1)

Base.length(S::AverageZero{Ulam})

Return the size of the Average Zero space

integral_pairing(ϕ::Observable{<:Ulam}, ρ::AbstractVector, ρ_w_error;
                 ρ_dual_weak_bound = …)

For Ulam the pairing simplifies: $ρ_N$ is piecewise constant and $ϕ.v[i] = N\,\int_{I_i} ϕ\,dx$ is N times the exact cell integral, so

\[\frac{1}{N}\sum_i ϕ.v[i]\,ρ_N[i] \;=\; \sum_i ρ_N[i]\,\int_{I_i} ϕ\,dx \;=\; \int_0^1 ϕ\,ρ_N\,dx.\]

There is no projection-error term $\|ϕ - ϕ_N\|_w\,\|ρ_N\|_{w^*}$ to add here — the discrete pairing already equals $\int ϕ\,ρ_N$ exactly (modulo floating-point rounding, which interval lifting of ρ covers). The only error contribution is from $ρ \neq ρ_N$: $|\int ϕ\,(ρ - ρ_N)\,dx| \leq \|ϕ\|_{L^∞}\,\|ρ - ρ_N\|_{L^1}$.

The ρ_dual_weak_bound kwarg is accepted for API symmetry with the Fourier method but is unused for Ulam.

nonzero_on(B::Ulam, (a, b))

Returns the indices of the elements of the Ulam basis that intersect with the interval y We do not assume an order of a and b; this should not matter unless the preimages are computed with very low accuracy. We assume, though, that y comes from the (possibly inexact) numerical approximation of an interval in $[0,1]$, i.e., we restrict to $y \cap [0,1]$

relative_measure((a,b)::Tuple{<:Interval,<:Interval}, (c,d)::Tuple{<:Interval,<:Interval})

Relative measure of the intersection of (a,b) wrt the whole interval (c,d) Assumes that a,b and c,d are sorted correctly

Hat

This type represents a Hat basis on $S^1$. It contains a vector with the midpoints of the hats.

Hat(n::Integer)

This constructs a Hat basis on $S^1$ on equispaced points

HatFunctionOnTorus

Hat function (on the torus)

This is a piecewise linear function such that: f(x) = 0 if x <= lo f(mi) = 1 f(x) = 0 if x >= ho

Evaluate a HatFunctionOnTorus correctly on an IntervalOnTorus

Assumption: this is only called on functions defined on our partition, so either mi==0, hi==0, or the three values are in increasing order

IntervalOnTorus

A separate type for intervals on the torus (mod 1) to "remind" us of the quotient

The interval is normalized in the constructor: the caller may assume that

• 0 <= inf(i) < 1
• sup(i) < inf(i) + 1 OR i==interval(0,1)

Base.getindex(B::Hat, i::Int)

Make so that B[j] returns a HatFunctionOnTorus with the j-th basis element

Given a preimage $y$ of a point $x$, this iterator returns $\phi_j(y)/T'(y)$

Base.length(B::Hat{T})

Return the size of the Hat basis

nonzero_on(B::Hat, dual_element)

Return the range of indices of the elements of the basis whose support intersects with the given dual element (i.e., a pair (y, absT')). The range may end with length(B)+1; this must be interpreted "mod length(B)": it means that it intersects with the hat function peaked in 0 as well (think for instance y = 0.9999).

Hat function (on the reals)

This is a piecewise linear function such that: f(x) = 0 if x <= lo f(mi) = 1 f(hi)

Evaluate a HatFunction (on the real line).

Must work correctly when S is an interval.

Base.getindex(B::HatNP, i::Int)

makes so that B[j] returns a HatFunction with the j-th basis element

Given a preimage y of a point x, this iterator returns \phi_j(y)/T'(y)

Base.length(B::HatNP)

Return the size of the HatNP basis

Return the range of indices of the elements of the basis whose support intersects with the given dual element (i.e., a pair (y, absT')). The range may end with length(B)+1; this must be interpreted "mod length(B)": it means that it intersects with the hat function peaked in 0 as well (think for instance y = 0.9999).

Equispaced partition of size n of [0,1]

Return (in an iterator) the pairs (i, (x, |T'(x)|)) where x is a preimage of p[i], which describe the "dual" L* evaluation(p[i])

Given a preimage y of a point x, this iterator returns $\phi_j(y)/T'(y)$

Return the size of the C2 basisBase.length(S::AverageZero) = length(S.basis)-1

Return the range of indices of the elements of the basis whose support intersects with the given dual element (i.e., a pair (y, absT', derder)).

Base.getindex(B::Chebyshev, i::Int)

Make so that B[j] returns a HatFunctionOnTorus with the j-th basis element

Return the size of the Chebyshev basis

_eval_T

Eval the Chebyshev polynomial up to degree n on an array of points in [-1, 1].

Not satisfactory, the intervals explode

bound_linalg_norm_L1_from_weak(B::Chebyshev)

Returns A such that ||ĉ||_{ℓ¹} ≤ A·||f||_{C1}.

Chebyshev coefficients satisfy |cⱼ| ≤ 2·||f||∞ for j ≥ 1 and |c₀| ≤ ||f||∞, so ||ĉ||{ℓ¹} ≤ (2n-1)·||f||∞ ≤ (2n-1)·||f||_{C1} where n = length(B).

bound_linalg_norm_L∞_from_weak(B::Chebyshev)

Returns A such that ||ĉ||_{ℓ∞} ≤ A·||f||_{C1}.

maxj |cⱼ| ≤ 2·||f||∞ ≤ 2·||f||_{C1}.

bound_weak_norm_from_linalg_norm(B::Chebyshev)

Returns (W₁, W₂) such that ||f||_{C1} ≤ W₁·||ĉ||_{ℓ¹} + W₂·||ĉ||_{ℓ∞}.

Since |Tⱼ(x)| ≤ 1: ||f||∞ ≤ ||ĉ||{ℓ¹}. For the derivative: ||f'||∞ ≤ Σ|cⱼ|·2j² ≤ 2(n-1)²·||ĉ||{ℓ¹} where n = degree.

certify_spectral_gap(B::Chebyshev, Q::DiscretizedOperator;
                      samples=256, radius_factor=0.5)

Run the full spectral gap certification pipeline for a Chebyshev discretized operator.

Uses BallArithmetic's CertifScripts to:

1. Compute Schur decomposition with rigorous error bounds
2. Find the spectral gap (distance from eigenvalue 1 to next largest eigenvalue)
3. Certify the resolvent on a circle separating eigenvalue 1 from the rest
4. Compute projector error bound δP for `restricttoaveragezero`

Returns a named tuple with fields:

• spectral_gap, second_largest, ρ (contour radius)
• M_inf (certified resolvent bound on contour)
• projector_error (δ_P for projector inflation)
• certification, schur_data, eigenvalue_index

eval_Clenshaw_BackwardFirst

Eval a polynomial in Chebyshev basis, ClenshawBackward, using ball arithmetic Following Viviane Ledoux, Guillaume Moroz "Evaluation of Chebyshev polynomials on intervals andapplication to root finding"

Return the norm of the matrix the NormCacher is working on.

restrict_to_average_zero(B::Chebyshev, BM::BallMatrix, f; certification=nothing)

Restrict BM to the average-zero subspace U⁰ using a certified spectral projector.

For Chebyshev, integral_covector(B) ≠ [1,0,...,0], so simple submatrix extraction does not work. Instead, we compute the Riesz projector P₁ onto the eigenspace of eigenvalue 1 (the invariant measure direction) via Schur decomposition, then return (I - P₁) * BM.

If certification is provided (from certify_spectral_gap), the projector radii are inflated by the certified projector error δ_P to account for the Schur factorization error ‖ZTZ* - A‖₂.

weak_by_strong_and_aux_bound(B::Chebyshev)

Returns (S₁, S₂) such that ||f||_{C1} ≤ S₁·||f||_{W^{k,1}} + S₂·||f||_{L1}.

For k ≥ 2: Sobolev embedding gives ||f||∞ ≤ ||f||{L1} + ||f'||{L1} and ||f'||∞ ≤ ||f'||{L1} + ||f''||{L1}, so ||f||{C1} ≤ 2·||f||{W^{k,1}}.

For k = 1: uses Markov inequality ||p'||∞ ≤ 2n²·||p||∞ for polynomials of degree n on [0,1], giving ||f||{C1} ≤ (1 + 2n²)·||f||{W^{1,1}}.

p1 * p2  for two `ProjectedFunction{<:Fourier}` on the same basis

Multiply two Fourier projections by treating each operand as a trigonometric polynomial (the discrete coefficient vector) plus an opaque L²-error bound. All bounds are derived from the finite coefficient vectors alone — provenance (the original W^{k,1} seminorm, function class, …) is not used. This keeps multiplication composable: the result is again a trigonometric polynomial with an L² error bound, ready to be passed to * or + again.

The arithmetic:

1. Convolve the two N-mode coefficient vectors → 2N−1 modes.
2. Truncate the central N modes back into the basis layout.
3. Parseval gives $\|\text{discarded modes}\|_{\ell^2}$ directly from the convolution output.

The bounds:

• $\|fg\|_{L^2} \leq \|\hat f\|_{\ell^1} \cdot \|g\|_{L^2}$ (Young's inequality, with $\|\hat f\|_{\ell^1}$ the surrogate for $\|f\|_{L^\infty}$).
• ``\|fg - πN(φf φg)\|{L^2} \leq \|\hat f\|{\ell^1}\,p2.\text{proj_error}
  • \|\hat g\|{\ell^1}\,p1.\text{proj_error} + \|\text{discarded modes}\|_{\ell^2}``.

Each input's proj_error is the L² distance from the original function to its trigonometric-polynomial approximant; the multiplication treats that distance as the only thing it knows.

Make so that B[j] returns the basis function of coordinate j

integral_pairing(ϕ::Observable{<:Fourier}, ρ, ρ_w_error;
                 ρ_dual_weak_bound = weak_dual_norm_bound(ϕ.B, ρ))

For Fourier bases (weak L²), the pairing $\int_0^1 ϕ_N(x)\,ρ_N(x)\,dx$ equals $\sum_n \hat ϕ_n\,\overline{\hat ρ_n}$. The result is real-valued for real signals; we return the real part of the sum.

Return the size of the Fourier basis

Return the weak-strong norm bound when restricted on the finite dimensional subspace. For Aη, by Cauchy-Schwarz: $||v||_{Aη} = \sum |ĉ_k| e^{2πη|k|} \leq (\sum e^{4πη|k|})^{1/2} \cdot ||ĉ||_{ℓ²}$ Geometric series: $\sum_{|k|\leq N} e^{4πη|k|} = 1 + 2(e^{4πη} - e^{4πη(N+1)})/(1 - e^{4πη})$

For W^{k,1}: Bernstein + Cauchy-Schwarz on [0,1]: $||v^{(j)}||_{L¹} ≤ ||v^{(j)}||_{L²} ≤ (2πN)^j ||v||_{L²}$ Sum: $M₁n = \sum_{j=0}^{k} (2πN)^j$