NonlinearIntegrators

Continuous Galerkin Variational Integrator.

• b: weights of the quadrature rule
• c: nodes of the quadrature rule
• x: nodes of the basis
• m: mass matrix
• a: derivative matrix
• r₀: reconstruction coefficients at the beginning of the interval
• r₁: reconstruction coefficients at the end of the interval

bias_grid(lo, hi, n, ::Type{T}) -> Vector{T}

Uniform grid of n + 1 bias values from lo to hi in precision T.

The coordinates are generated from an integer-indexed range in Float64 and then cast to T, so a large n cannot overflow the step to zero in reduced precision (the Float16 trap where T(n) overflows to Inf and (hi - lo)/n evaluates to zero, previously throwing ArgumentError: range step cannot be zero). Only the grid coordinates touch Float64; the seed's dictionary and solve run entirely in T.

initial_params!(int, ::OGA1d_Legacy)

Legacy OGA initial guess for NonLinear_OneLayer_GML, kept as a selectable alternative to the default OGA1d for comparison. This is the pre-refactor algorithm: the dictionary and the greedy least-squares fit are assembled in Float64 (a "double-precision island"), the output weights are obtained from the normal equations Gk \ rhs, and the result is rounded into the working-precision cache. See the "Orthogonal Greedy Algorithm" section of the documentation for why this was replaced by a working-precision QR fit. Select it with NonLinear_OneLayer_GML(...; initial_guess_method = OGA1d_Legacy()).

initial_trajectory!(sol, history, params, int, initial_trajectory)

Initial trajectory for the NonLinear_OneLayer_Lux integrator.

oga_norm_floor(::Type{T}, ref) -> T

Amplitude-scale floor below which a dictionary-atom norm is treated as numerically zero (so normalization is skipped rather than dividing by noise). Scales with sqrt(eps(T)) * ref, i.e. relative to the largest atom ref and to the working precision: a norm smaller than this cannot be inverted reliably in T.

Replaces the hard-coded absolute 1e-12 guard, which sat below eps(Float32) and so never fired in reduced precision.

oga_tikhonov(G; C = 100) -> eltype(G)

Scale- and precision-relative Tikhonov floor C * eps(T) * tr(G)/n for stabilizing a Gram / normal-equations solve G \ rhs at working precision T = eltype(G). Relative to the mean diagonal tr(G)/n, so the effective condition-number cap tracks the precision; C is a modest safety factor.

Replaces the hard-coded absolute 1e-12 / 1e-14 ridges, which round away entirely below eps(Float32). Provided as a fallback for a regularized normal-equations variant; the OGA fit itself uses weighted_lstsq, which avoids forming G.

weighted_lstsq(Φ, w, y; C = 100) -> Vector

Solve the (ridged) quadrature-weighted least-squares fit for the output weights x,

minₓ  Σⱼ wⱼ (Σᵢ xᵢ Φ[i,j] − yⱼ)²  +  λ‖x‖² ,

where each row of Φ is a dictionary atom sampled at the quadrature nodes, w are the (positive) quadrature weights, and y is the fit target at those nodes.

Solved by QR on the √w-scaled design matrix rather than the normal equations Φ diag(w) Φᵀ, so accuracy is governed by κ(Φ) instead of κ(Φ)² — the difference that lets the fit run in reduced precision without a rank-deficient Gram matrix.

The plain QR solve is used whenever it is finite (so the Float64/Float32 atom choice is byte-for-byte the old Gram solution and the greedy residual is unperturbed). Only if it returns a non-finite result — a genuinely rank-deficient design matrix, i.e. the Float16 case — is the fit retried with a √λ·I augmentation, where the ridge λ = C · eps(T) · tr(ÂᵀÂ)/natoms is the precision-scaled Tikhonov floor (see oga_tikhonov) that keeps the solution bounded and the seed finite.