CohomologicalEquations

CohomologicalEquations

Solve the cohomological equations that arise in the parametrisation method for computing Spectral Submanifolds (SSMs) of high-dimensional dynamical systems.

Problem statement

Given a full‑order model of order ORD in FOM degrees of freedom, the parametrisation method seeks a change of coordinates

x(t) = W(z(t)),          z ∈ ℂᴺᵛᵃʳ

such that the reduced dynamics ż = R(z) is simpler. Expanding both W and R as dense polynomials in the NVAR = ROM + N_EXT reduced variables (ROM master‑mode amplitudes plus N_EXT external forcing amplitudes) and matching coefficients monomial by monomial yields the cohomological equations: one linear system per multi‑index α.

System structure for multi‑index α

Let s = ⟨λ, α⟩ be the superharmonic frequency and let nR = |ℛ(α)| be the number of master modes resonant at α. The cohomological linear system is

┌              ┐ ┌         ┐   ┌           ┐
│  L(s)  C(s)  │ │  W[α]   │ = │  RHS_inv  │   FOM rows  (invariance)
│  L̂(s)  Ĉ(s)  │ │  R_res  │   │  RHS_ort  │   nR  rows  (orthogonality)
└              ┘ └         ┘   └           ┘

where:

• L(s) (FOM × FOM) is the parametrisation operator,
• C(s) (FOM × nR) acts on the unknown resonant reduced-dynamics coefficients,
• L̂(s) (nR × FOM) is the orthogonality row operator,
• Ĉ(s) (nR × nR) is the orthogonality joint operator,
• W[α] ∈ ℂᶠᵒᵐ is the (zeroth-order) parametrisation coefficient,
• R_res ∈ ℂⁿᴿ are the reduced-dynamics coefficients for resonant modes.

Non-resonant master modes are trivially zero and are excluded from the system. External forcing modes are known and appear only on the right-hand side.

Module contents

CohomologicalContext{T, ORD, ORDP1, NVAR, ROM, FOM}

A single flat struct that bundles all precomputed data required to solve the cohomological equations for every monomial in the parametrisation method.

Using one flat struct rather than nested containers eliminates the name collision between the invariance‑equation operators and the orthogonality‑condition operators (both previously called C_coeffs/E_coeffs) and makes the data provenance explicit at every call site.

Type parameters

Fields

Full‑order linear operators

• linear_terms :: NTuple{ORDP1, Matrix{T}} — coefficient matrices of the linear part of the full‑order model, (B₀, B₁, …, B_ORD), each of size FOM × FOM.

• generalised_eigenmodes :: Matrix{T} — size FOM × NVAR. Right generalised eigenvectors; columns 1:ROM are the master modes, columns ROM+1:NVAR are the external forcing modes.

• lambda_diag :: Vector{T} — length NVAR. Diagonal entries λᵢ of the Jordan matrix, used to form the superharmonic frequency s = ⟨λ, α⟩ for multi‑index α. These are read directly from the reduced‑dynamics polynomial R at the linear monomials e_i.

Invariance‑equation operators (polynomial coefficients in s)

These arrays are produced by InvarianceEquation.precompute_column_polynomials.

• invariance_C_coeffs :: Vector{Matrix{T}} — length ROM. invariance_C_coeffs[r][:, j] is the degree‑(j-1) coefficient of the column operator C_r(s) acting on the reduced‑dynamics unknown for master mode r; each matrix has size FOM × ORD.

• invariance_E_coeffs :: Vector{Matrix{T}} — length N_EXT. invariance_E_coeffs[e][:, j] is the degree‑(j-1) coefficient of the external‑forcing operator E_e(s) for external variable e; each matrix has size FOM × ORD.

Orthogonality‑condition operators (polynomial coefficients in s)

These arrays are produced by MasterModeOrthogonality.precompute_orthogonality_operator_coefficients and MasterModeOrthogonality.precompute_orthogonality_column_polynomials.

• orthogonality_J_coeffs :: Vector{Matrix{T}} — length ROM. orthogonality_J_coeffs[r][j, :] is the degree‑(j-1) coefficient of the left row operator L_r(s) for master mode r; each matrix has size ORD × FOM.

• orthogonality_C_coeffs :: Vector{Matrix{T}} — length ROM. orthogonality_C_coeffs[r][j, :] is the degree‑(j-1) coefficient of the joint column operator Ĉ_r(s) acting on the master‑mode unknowns; each matrix has size (ORD-1) × ROM.

• orthogonality_E_coeffs :: Vector{Matrix{T}} — length ROM. orthogonality_E_coeffs[r][j, :] is the degree‑(j-1) coefficient of the joint column operator Ê_r(s) acting on the external variables; each matrix has size (ORD-1) × N_EXT.

Resonance

• resonance_set :: ResonanceSet — look‑up table indicating which master modes are resonant with each monomial.

Precomputed bookkeeping

• linear_monomial_skip_set :: Set{Int} — indices of all unit-vector monomials (e_1, …, e_NVAR); these are initialised before the main loop and skipped.

solve_cohomological_equations!(W, R, ctx, model) -> nothing

Solve the cohomological equations for all monomials in the multiindex set of W and R, processing them in causal order (ascending total degree so that lower‑order coefficients are available when higher‑order ones are solved).

All unit‑vector monomials eᵣ for r = 1 … NVAR (both master modes and external forcing modes) are assumed to have been initialised to the spectral data beforehand and are skipped.

Arguments

• W :: Parametrisation{ORD, NVAR, T} — parametrisation (updated in‑place).
• R :: ReducedDynamics{ROM, NVAR, T} — reduced dynamics (updated in‑place).
• ctx :: CohomologicalContext — all precomputed data and the resonance set.
• model :: NDOrderModel — full‑order model; passed through to solve_single_monomial! for nonlinear term evaluation.
• ml_cache :: MultilinearTermsCache — precomputed factorisation cache built once before the loop.

solve_cohomological_problem(
	model, mset, master_eigenvalues,
	master_modes, left_eigenmodes, resonance_set;
	initial_W = nothing, initial_R = nothing
) -> (W, R)

High‑level driver that assembles a CohomologicalContext from raw spectral data and solves the full set of cohomological equations.

External eigenvalues are read directly from model.external_system.eigenvalues (or treated as absent when model.external_system === nothing). The linear‑operator tuple is read from model.linear_terms.

Steps

1. Extract external eigenvalues from model.external_system and build the Jordan matrix Λ.
2. Create (or reuse) the Parametrisation W and ReducedDynamics R objects and initialise the master‑mode linear monomials from master_modes.
3. Solve the linear cohomological equations for each external forcing direction: L(s_ext) · W_ext = compute_multilinear_terms(model, e_ext, W). This yields the particular solution at the forcing frequency and populates the external columns of generalised_right_eigenmodes.
4. Build the N_EXT × L external‑dynamics matrix (non‑zero only at the linear forcing monomials e_{ROM+e}).
5. Precompute the invariance‑equation operator columns (invariance_C_coeffs, invariance_E_coeffs) via InvarianceEquation.precompute_column_polynomials.
6. Precompute the orthogonality‑condition operators (orthogonality_J_coeffs, orthogonality_C_coeffs, orthogonality_E_coeffs) via MasterModeOrthogonality.precompute_orthogonality_operator_coefficients and MasterModeOrthogonality.precompute_orthogonality_column_polynomials.
7. Assemble a CohomologicalContext and call solve_cohomological_equations!.

Arguments

• model :: NDOrderModel — full‑order model. model.linear_terms provides (B₀, B₁, …, B_ORD); model.external_system (may be nothing) carries the external eigenvalues.
• mset :: MultiindexSet{NVAR} — multiindex set over all NVAR reduced variables (NVAR = ROM + N_EXT).
• master_eigenvalues :: SVector{ROM, ComplexF64} — eigenvalues of the master modes.
• master_modes :: Matrix{ComplexF64} — size FOM × ROM; right eigenvectors of the master modes. The external forcing directions are derived internally by solving the linear cohomological equations.
• left_eigenmodes :: AbstractMatrix{ComplexF64} — left eigenvectors of the master modes, size FOM × ROM (used in the orthogonality conditions).
• resonance_set :: ResonanceSet — precomputed resonance look‑up table.
• initial_W, initial_R — optionally supply already-initialised objects (their linear monomials must be set correctly).

Returns

(W, R) — the solved Parametrisation and ReducedDynamics.

solve_single_monomial!(W, R, idx, ctx, model) -> nothing

Solve the cohomological equations for the monomial with multiindex‑set position idx, updating the coefficients of W and R in‑place.

Algorithm

1. Compute the superharmonic frequency s = ⟨λ, α⟩ from the diagonal of ctx.lambda_diag.
2. Look up the resonance pattern resonance ∈ {true,false}^{ROM} from ctx.resonance_set.
3. Compute lower‑order coupling vectors ξ[j] (length FOM) via LowerOrderCouplings.compute_lower_order_couplings.
4. Compute the nonlinear model RHS in-place via MultilinearTerms.compute_multilinear_terms! into ctx.ml_result_buffer. No heap allocation occurs.
5. Retrieve the known external dynamics at this monomial from the last N_EXT rows of R.poly.coefficients[:, idx].
6. Assemble the stacked (FOM + nR) × (FOM + nR) linear system directly into ctx.system_matrix_buffer and ctx.rhs_buffer via the in-place variants InvarianceEquation.assemble_cohomological_matrix_and_rhs! and MasterModeOrthogonality.assemble_orthogonality_matrix_and_rhs!, then factor ctx.system_matrix_buffer in-place with lu!(view(...), check=false) and solve with ldiv!. The buffer is never reused after factorisation, so no copy is needed.
7. Store W[α] (zeroth time‑derivative) and the resonant reduced‑dynamics coefficients R_res.
8. Compute higher time‑derivative coefficients W^(j)[α] for j = 1 … ORD-1 via ParametrisationMethod.compute_higher_derivative_coefficients!.

Arguments

• W :: Parametrisation{ORD, NVAR, T} — parametrisation object (updated in‑place).
• R :: ReducedDynamics{ROM, NVAR, T} — reduced dynamics object (updated in‑place).
• idx :: Int — position of the target monomial in the shared multiindex set.
• ctx :: CohomologicalContext{T, ORD, ORDP1, NVAR, FOM} — all precomputed data.
• model :: NDOrderModel — full‑order model; provides nonlinear term evaluations.
• ml_cache :: MultilinearTermsCache — precomputed factorisation cache; built once before the solve loop.