Invariant solutions¶

Equilibria, travelling waves and periodic orbits found by Newton–Krylov–hookstep on the flow map of your own time stepper.

An invariant solution is a state that the flow maps to itself, up to a symmetry: an equilibrium \(v(u) = 0\), a travelling wave \(u(x - ct)\), a relative periodic orbit \(\tau_\ell\,\Phi_T(u) = u\). Conatus solves these equations with a Newton–Krylov method and a trust-region (hookstep) step. The Jacobian is never formed: each Krylov vector costs one linearized evaluation or time integration, exact for the discrete scheme.

The flow map¶

The searches use the time stepper of the model, here Kuramoto–Sivashinsky on a periodic domain of length 22:

import numpy as np

import conatus as cn
from conatus.expert import (NewtonKrylovHookstep, equilibrium_eigenvalues, find_equilibrium,
                            find_travelling_wave)

L = 22.0
m = cn.Model(cn.Box(x=cn.periodic(L)), name="Kuramoto-Sivashinsky")
m.scalar("u")
m.equations("dt(u) + lap(u) + lap(lap(u)) = -u*dx(u)")
sim = m.simulate(N=32, scheme="RK443", dt=0.01)        # the flow map used by the searches

newton = NewtonKrylovHookstep(tol=1e-12)    # accuracy required; algorithm controls defaulted
print(newton.explain())                     # budget, GMRES, trust radii, marked "(default)"

The tolerance of the Newton–Krylov solver is required (it decides the digits of the solution); its algorithm controls (iteration budget, Krylov tolerance and size, initial, minimum and maximum trust radius) have declared defaults, printed by newton.explain() and marked “(default)”: they change how fast the residual reaches the tolerance, not the answer, whose residual is returned.

Scope. find_equilibrium, find_travelling_wave and find_relative_periodic_orbit work on problems with one periodic direction (a real Fourier basis, no other axis): Kuramoto–Sivashinsky, 1D Swift–Hohenberg and similar. They accept the Simulation (or its sim.solver; a cn.Model with N= for equilibria and waves, whose search does not depend on the time scheme). Anything else is refused with this scope and the alternatives: steady states of large problems in symmetry subspaces (plane Couette, conatus.ecs.steady.AffineSteadyProblem), pipe travelling waves (conatus.ecs.pipe.PipeTravellingWave, experimental: Experimental features), 2D Kolmogorov flow (conatus.ecs.kolmogorov). A generalisation to every geometry is planned, not available.

The seed x0 of every finder is the state vector of the solver (sim.solver.X), the simulation itself (its current state), or the grid values of the unknown (sim.grid("u")), converted to the state vector; the conversion and the size of the round trip are in the result’s notes, and any other array is refused with the accepted shapes. The stability calls (equilibrium_eigenvalues, floquet_multipliers) take the simulation as the finders do. A Newton search converges only from a seed near the solution: most states of a chaotic run fall back to an equilibrium or fail (tutorial 5 keeps stored seeds). A result prints on one line (print(tw)); explain() has the evidence.

Equilibria¶

q = 2 * np.pi / L
sim.initial(u=lambda x: np.sin(q * x))                 # odd seed: u(-x) = -u(x)
E1 = find_equilibrium(sim, sim.solver.X.copy(), newton)   # the Simulation itself
print(E1.explain())
print(equilibrium_eigenvalues(sim.solver, E1, n_modes=4, dense_max=500, arnoldi_tol=1e-12))

An equilibrium of a translation-invariant system is not isolated (every translate is one), so find_equilibrium works in the subspace of odd states \(u(-x) = -u(x)\), where translations are broken; the seed must lie in it, which is checked. equilibrium_eigenvalues returns the leading eigenvalues of the Jacobian in the full space, translation mode included.

Travelling waves¶

sim.initial(u=lambda x: np.cos(q * x) + np.sin(2 * q * x))
TW1 = find_travelling_wave(sim.solver, sim.solver.X.copy(), 0.7, newton)
print(f"c = {TW1.params['c']:.5f}, residual {TW1.residual:.1e}")

The speed c is an unknown, fixed by a phase condition. The residual is \(\|v(U) + c\,\partial_x U\| / \|U\|\), evaluated in the full space. A solution whose speed term \(|c|\,\|\partial_x U\|\) is below the Newton accuracy (TW_SPEED_TOL times tol times the scale of the solve) is an equilibrium, not a travelling wave: it is returned as an Equilibrium (kind "equilibrium", params["c"] kept, the measured value in the notes). A Newton run that lands on \(u = 0\) is retracted as the trivial state (TRIVIAL_TOL); the boundary-value NewtonSolver flags the same case with trivial=True and a note in explain() (it is not refused there: \(u = 0\) may be the wanted answer).

Periodic and pre-periodic orbits¶

find_relative_periodic_orbit(solver, X0, T0, ell0, newton, steps=..., reflect=False) solves \(\tau_\ell\,\Phi_T(u) = u\) for \((u, T, \ell)\) with time and translation phase conditions; reflect=True gives pre-periodic orbits \(R\,\tau_\ell\,\Phi_T(u) = u\). The number of time steps over one period is part of the problem. floquet_multipliers returns the multipliers of the orbit and identifies its marginal directions. Seeds come from near-recurrences of a long run; tutorial 5. Invariant solutions: Kuramoto–Sivashinsky, L = 22 carries one to the published Floquet exponents.

Thresholds and edge states¶

The boundary between initial states that decay and those that become turbulent is located by bisection on an amplitude. Conatus returns the threshold as a bracket and never counts an undecided run on either side:

from conatus.expert import ThresholdBisection

def fate(a):
    """Declared verdict for amplitude a: x' = -x (x - 1)(x - 2), x(0) = a, t in [0, 40]."""
    x = a
    for _ in range(4000):
        x += 0.01 * (-x * (x - 1) * (x - 2))
    return -1 if x < 0.5 else (+1 if x > 1.5 else 0)

tb = ThresholdBisection(fate, rtol=1e-4)
lo, hi = tb.run(0.5, 1.6)                    # the threshold as a bracket ]lo, hi]
print(tb.explain())

The verdict function is yours (a time integration and declared thresholds). EdgeTracker repeats the bisection along the trajectory to follow the edge state itself. For 3D shear flows, NLOP.from_solver (module conatus.adjoint.nlop) computes the nonlinear optimal perturbation that seeds such a bisection.

Validation¶

A solution is accepted on its residual, reported with the iteration history (explain()). Its stability is a separate computation (eigenvalues or Floquet multipliers in the full space). A solution found in a symmetry subspace may be unstable outside it (see Symmetry subspaces).