Skip to content

Reference

Analysis

The quantifiers that consume any System. Prose-first treatments live in the Analysis section.

Orbit diagrams

orbit_diagram 

orbit_diagram(
    sys: Any,
    param: str,
    values: Any,
    *,
    n: int = 200,
    transient: int = 500,
    carry_state: bool = True,
    components: int | str | tuple = 0,
    ic: Any | None = None,
) -> OrbitDiagram

Sweep a parameter and record the asymptotic orbit at each value.

Works on anything discrete: a :class:~tsdynamics.families.DiscreteMap directly, or a flow wrapped in a :class:~tsdynamics.derived.PoincareMap / :class:~tsdynamics.derived.StroboscopicMap — in which case this is the bifurcation diagram of the flow. ODE parameter changes reuse the compiled module (control parameters), so flow sweeps stay cheap; DDE sweeps recompile per value (their structure depends on all parameters).

PARAMETER DESCRIPTION
sys

The system to sweep. Never mutated — each value gets a fresh with_params copy.

TYPE: System(discrete)

param

Parameter name to sweep.

TYPE: str

values

Parameter values, in sweep order.

TYPE: iterable of float

n

Points recorded per parameter value.

TYPE: int DEFAULT: 200

transient

Steps discarded before recording, at every value.

TYPE: int DEFAULT: 500

carry_state

Start each value from the previous value's final state (follows the attractor branch; the classic way to draw clean diagrams). When False, every value starts from ic / the system default.

TYPE: bool DEFAULT: True

components

Which state components to record (names allowed when the system declares variables).

TYPE: int, str, or tuple DEFAULT: 0

ic

Initial state for the first value (and every value when carry_state=False).

TYPE: array - like DEFAULT: None

Examples:

>>> od = orbit_diagram(Logistic(), "r", np.linspace(2.5, 4.0, 600), n=120)
>>> x, y = od.flat()
>>> # bifurcation diagram of a flow:
>>> od = orbit_diagram(PoincareMap(Rossler(), (1, 0.0)), "c", np.linspace(2, 6, 80))
Source code in src/tsdynamics/analysis/orbits/orbit_diagram.py 
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
def orbit_diagram(
    sys: Any,
    param: str,
    values: Any,
    *,
    n: int = 200,
    transient: int = 500,
    carry_state: bool = True,
    components: int | str | tuple = 0,
    ic: Any | None = None,
) -> OrbitDiagram:
    """
    Sweep a parameter and record the asymptotic orbit at each value.

    Works on anything discrete: a :class:`~tsdynamics.families.DiscreteMap`
    directly, or a flow wrapped in a
    :class:`~tsdynamics.derived.PoincareMap` /
    :class:`~tsdynamics.derived.StroboscopicMap` — in which case this *is*
    the bifurcation diagram of the flow.  ODE parameter changes reuse the
    compiled module (control parameters), so flow sweeps stay cheap; DDE
    sweeps recompile per value (their structure depends on all parameters).

    Parameters
    ----------
    sys : System (discrete)
        The system to sweep.  Never mutated — each value gets a fresh
        ``with_params`` copy.
    param : str
        Parameter name to sweep.
    values : iterable of float
        Parameter values, in sweep order.
    n : int
        Points recorded per parameter value.
    transient : int
        Steps discarded before recording, at every value.
    carry_state : bool
        Start each value from the previous value's final state (follows the
        attractor branch; the classic way to draw clean diagrams).  When
        False, every value starts from ``ic`` / the system default.
    components : int, str, or tuple
        Which state components to record (names allowed when the system
        declares ``variables``).
    ic : array-like, optional
        Initial state for the first value (and every value when
        ``carry_state=False``).

    Examples
    --------
    >>> od = orbit_diagram(Logistic(), "r", np.linspace(2.5, 4.0, 600), n=120)
    >>> x, y = od.flat()
    >>> # bifurcation diagram of a flow:
    >>> od = orbit_diagram(PoincareMap(Rossler(), (1, 0.0)), "c", np.linspace(2, 6, 80))
    """
    if not sys.is_discrete:
        raise TypeError(
            "orbit_diagram needs a discrete-time view: a DiscreteMap, or a flow wrapped "
            "in PoincareMap / StroboscopicMap."
        )

    comp = (components,) if isinstance(components, int | str) else tuple(components)
    # Resolve names via the *instance* (not ``type(sys)``): a derived wrapper
    # exposes ``variables`` as a property, so ``type(sys).variables`` returns the
    # descriptor object (truthy) and short-circuits — breaking named components
    # over a PoincareMap/StroboscopicMap.  Instance lookup returns the ClassVar
    # for families and the resolved names for wrappers alike.
    names = getattr(sys, "variables", None)
    idx: list[int] = []
    for c in comp:
        if isinstance(c, str):
            if names is None:
                raise ValueError("named components need the system to declare `variables`")
            idx.append(names.index(c))
        else:
            idx.append(int(c))

    import warnings

    values_arr = np.asarray(list(values), dtype=float)
    points: list[np.ndarray] = []
    state: np.ndarray | None = None

    for v in values_arr:
        current = sys.with_params(**{param: v})
        start = state if (carry_state and state is not None) else ic
        try:
            current.reinit(start)
            for _ in range(transient):
                current.step()
            rec = np.empty((n, len(idx)))
            for i in range(n):
                u = current.step()
                rec[i] = u[idx]
        except RuntimeError as exc:
            # One divergent value must not discard the whole sweep: record an
            # empty point set and restart the next value from `ic`.
            warnings.warn(
                f"orbit_diagram: {param}={v:g} diverged ({exc}); recording an "
                f"empty set for this value.",
                RuntimeWarning,
                stacklevel=2,
            )
            points.append(np.empty((0, len(idx))))
            state = None
            continue
        points.append(rec)
        if carry_state:
            state = current.state()

    meta = {
        "system": type(sys).__name__,
        "param": param,
        "n": n,
        "transient": transient,
        "carry_state": carry_state,
        "components": tuple(idx),
    }
    return OrbitDiagram(
        param=param, values=values_arr, points=points, components=tuple(idx), meta=meta
    )

OrbitDiagram dataclass 

OrbitDiagram(
    param: str,
    values: ndarray,
    points: list[ndarray],
    components: tuple[int, ...],
    meta: dict = dict(),
)

Result of :func:orbit_diagram.

Iterate to get (value, points) pairs, or use :meth:flat for the scatter-ready arrays.

flat 

flat(component: int = 0) -> tuple[ndarray, ndarray]

Flatten to scatter-plot arrays (x, y).

x repeats each parameter value once per recorded point; y is the chosen recorded component.

Source code in src/tsdynamics/analysis/orbits/orbit_diagram.py 
34
35
36
37
38
39
40
41
42
43
44
45
def flat(self, component: int = 0) -> tuple[np.ndarray, np.ndarray]:
    """
    Flatten to scatter-plot arrays ``(x, y)``.

    ``x`` repeats each parameter value once per recorded point; ``y`` is
    the chosen recorded component.
    """
    x = np.concatenate(
        [np.full(p.shape[0], v) for v, p in zip(self.values, self.points, strict=True)]
    )
    y = np.concatenate([p[:, component] for p in self.points])
    return x, y

Poincaré sections

poincare_section 

poincare_section(
    sys_or_traj: Any,
    plane: tuple,
    *,
    direction: int = +1,
    steps: int = 1000,
    transient: int = 0,
    dt: float = 0.01,
    max_time: float = 10000.0,
) -> Trajectory

Poincaré surface of section.

Two input modes:

  • System → wraps it in a :class:~tsdynamics.derived.PoincareMap and collects steps root-refined crossings (accurate).
  • Trajectory → finds the plane crossings between consecutive samples by linear interpolation (pure data path; accuracy limited by the trajectory's sampling interval).
PARAMETER DESCRIPTION
sys_or_traj

TYPE: System or Trajectory

plane

(i, c) for the section y_i = c, or (normal, offset).

TYPE: tuple

direction

Crossing direction filter.

TYPE: (+1, -1, 0) DEFAULT: +1

steps

System mode only — see :class:~tsdynamics.derived.PoincareMap.

TYPE: int DEFAULT: 1000

transient

System mode only — see :class:~tsdynamics.derived.PoincareMap.

TYPE: int DEFAULT: 1000

dt

System mode only — see :class:~tsdynamics.derived.PoincareMap.

TYPE: int DEFAULT: 1000

max_time

System mode only — see :class:~tsdynamics.derived.PoincareMap.

TYPE: int DEFAULT: 1000

RETURNS DESCRIPTION
Trajectory

t = crossing times, y = full-dimensional crossing states.

Examples:

>>> section = poincare_section(Rossler(), plane=(0, 0.0), steps=500)
>>> section = poincare_section(traj, plane=(2, 25.0))     # from data
Source code in src/tsdynamics/analysis/orbits/poincare.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def poincare_section(
    sys_or_traj: Any,
    plane: tuple,
    *,
    direction: int = +1,
    steps: int = 1000,
    transient: int = 0,
    dt: float = 0.01,
    max_time: float = 1e4,
) -> Trajectory:
    """
    Poincaré surface of section.

    Two input modes:

    - **System** → wraps it in a :class:`~tsdynamics.derived.PoincareMap`
      and collects ``steps`` root-refined crossings (accurate).
    - **Trajectory** → finds the plane crossings between consecutive samples
      by linear interpolation (pure data path; accuracy limited by the
      trajectory's sampling interval).

    Parameters
    ----------
    sys_or_traj : System or Trajectory
    plane : tuple
        ``(i, c)`` for the section ``y_i = c``, or ``(normal, offset)``.
    direction : {+1, -1, 0}
        Crossing direction filter.
    steps, transient, dt, max_time
        System mode only — see :class:`~tsdynamics.derived.PoincareMap`.

    Returns
    -------
    Trajectory
        ``t`` = crossing times, ``y`` = full-dimensional crossing states.

    Examples
    --------
    >>> section = poincare_section(Rossler(), plane=(0, 0.0), steps=500)
    >>> section = poincare_section(traj, plane=(2, 25.0))     # from data
    """
    if isinstance(sys_or_traj, Trajectory):
        return _section_from_data(sys_or_traj, plane, direction)
    pmap = PoincareMap(sys_or_traj, plane, direction=direction, dt=dt, max_time=max_time)
    return pmap.trajectory(steps, transient=transient)

Lyapunov quantifiers

lyapunov_spectrum 

lyapunov_spectrum(sys: Any, **kwargs: Any) -> ndarray

Lyapunov spectrum of any system — uniform entry point.

Dispatches to the family implementation (QR tangent dynamics for maps, compiled variational equations for ODEs, jitcdde_lyap for DDEs). Keyword arguments are forwarded (steps= for maps; final_time=, dt=, burn_in=, ... for flows).

Source code in src/tsdynamics/analysis/lyapunov/__init__.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def lyapunov_spectrum(sys: Any, **kwargs: Any) -> np.ndarray:
    """
    Lyapunov spectrum of any system — uniform entry point.

    Dispatches to the family implementation (QR tangent dynamics for maps,
    compiled variational equations for ODEs, ``jitcdde_lyap`` for DDEs).
    Keyword arguments are forwarded (``steps=`` for maps; ``final_time=``,
    ``dt=``, ``burn_in=``, ... for flows).
    """
    method = getattr(sys, "lyapunov_spectrum", None)
    if method is None:
        raise TypeError(
            f"{type(sys).__name__} has no lyapunov_spectrum implementation. "
            f"For derived wrappers, compute the spectrum on the underlying system."
        )
    return method(**kwargs)

max_lyapunov 

max_lyapunov(
    sys: Any,
    *,
    d0: float = 1e-09,
    n_rescale: int = 400,
    steps_per: int = 5,
    dt: float | None = None,
    transient: int = 500,
    ic: Any | None = None,
    seed: int | None = None,
) -> float

Maximal Lyapunov exponent by two-trajectory rescaling (Benettin et al. 1976).

Runs a reference and a perturbed copy of the system in lockstep through the :class:~tsdynamics.families.System protocol — no Jacobian needed, so it works for any ODE or map (including ones with non-smooth right-hand sides). Not available for DDEs (their state cannot be set_state-ed); use DelaySystem.lyapunov_spectrum instead.

PARAMETER DESCRIPTION
sys

ODE or map.

TYPE: System

d0

Perturbation size restored at every rescaling.

TYPE: float DEFAULT: 1e-09

n_rescale

Number of rescaling cycles (more → better averaging).

TYPE: int DEFAULT: 400

steps_per

Protocol steps between rescalings.

TYPE: int DEFAULT: 5

dt

Step size for continuous systems (default: the system's step default).

TYPE: float DEFAULT: None

transient

Protocol steps discarded before measuring.

TYPE: int DEFAULT: 500

ic

Initial condition for the reference trajectory.

TYPE: array - like DEFAULT: None

seed

Seed for the random perturbation direction.

TYPE: int DEFAULT: None

RETURNS DESCRIPTION
float

Estimated maximal exponent (per unit time / per iteration).

Examples:

>>> max_lyapunov(Lorenz(ic=[1.0, 1.0, 1.0]), dt=0.05)   # ≈ 0.91
Source code in src/tsdynamics/analysis/lyapunov/__init__.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def max_lyapunov(
    sys: Any,
    *,
    d0: float = 1e-9,
    n_rescale: int = 400,
    steps_per: int = 5,
    dt: float | None = None,
    transient: int = 500,
    ic: Any | None = None,
    seed: int | None = None,
) -> float:
    """
    Maximal Lyapunov exponent by two-trajectory rescaling (Benettin et al. 1976).

    Runs a reference and a perturbed copy of the system in lockstep through
    the :class:`~tsdynamics.families.System` protocol — no Jacobian needed, so it
    works for any ODE or map (including ones with non-smooth right-hand
    sides).  Not available for DDEs (their state cannot be ``set_state``-ed);
    use ``DelaySystem.lyapunov_spectrum`` instead.

    Parameters
    ----------
    sys : System
        ODE or map.
    d0 : float
        Perturbation size restored at every rescaling.
    n_rescale : int
        Number of rescaling cycles (more → better averaging).
    steps_per : int
        Protocol steps between rescalings.
    dt : float, optional
        Step size for continuous systems (default: the system's step default).
    transient : int
        Protocol steps discarded before measuring.
    ic : array-like, optional
        Initial condition for the reference trajectory.
    seed : int, optional
        Seed for the random perturbation direction.

    Returns
    -------
    float
        Estimated maximal exponent (per unit time / per iteration).

    Examples
    --------
    >>> max_lyapunov(Lorenz(ic=[1.0, 1.0, 1.0]), dt=0.05)   # ≈ 0.91
    """
    if isinstance(sys, DelaySystem):
        raise NotImplementedError(
            "max_lyapunov needs set_state, which delay systems cannot support — "
            "use DelaySystem.lyapunov_spectrum (jitcdde_lyap) instead."
        )
    if sys.is_discrete and dt is not None:
        raise ValueError(
            "dt has no meaning for discrete maps — omit it (every step is one iteration)."
        )

    rng = np.random.default_rng(seed)
    ref = sys.copy()
    ref.reinit(ic)
    for _ in range(transient):
        ref.step(dt)

    pert = sys.copy()
    direction = rng.normal(size=sys.dim)
    direction *= d0 / np.linalg.norm(direction)
    pert.reinit(ref.state() + direction)

    log_sum = 0.0
    for _ in range(n_rescale):
        for _ in range(steps_per):
            ref.step(dt)
            pert.step(dt)
        delta = pert.state() - ref.state()
        d = float(np.linalg.norm(delta))
        if d == 0.0 or not np.isfinite(d):
            raise RuntimeError(
                "max_lyapunov: trajectories collapsed or diverged — "
                "try a larger d0 or smaller steps_per."
            )
        log_sum += np.log(d / d0)
        pert.set_state(ref.state() + (d0 / d) * delta)

    if sys.is_discrete:
        elapsed_per_cycle = float(steps_per)
    else:
        step_dt = dt if dt is not None else getattr(sys, "_default_step_dt", 0.01)
        elapsed_per_cycle = steps_per * float(step_dt)
    return log_sum / (n_rescale * elapsed_per_cycle)

kaplan_yorke_dimension 

kaplan_yorke_dimension(spectrum: Any) -> float

Kaplan–Yorke (Lyapunov) dimension from a Lyapunov spectrum.

D_KY = j + (λ₁ + ... + λ_j) / |λ_{j+1}| where j is the largest index with a non-negative cumulative sum (Kaplan & Yorke 1979).

PARAMETER DESCRIPTION
spectrum

Lyapunov exponents (any order; sorted descending internally).

TYPE: array - like

RETURNS DESCRIPTION
float

0.0 when every exponent is negative; len(spectrum) when the cumulative sum never turns negative (spectrum incomplete).

Examples:

>>> kaplan_yorke_dimension([0.906, 0.0, -14.57])   # Lorenz
2.062...
Source code in src/tsdynamics/analysis/lyapunov/__init__.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def kaplan_yorke_dimension(spectrum: Any) -> float:
    """
    Kaplan–Yorke (Lyapunov) dimension from a Lyapunov spectrum.

    ``D_KY = j + (λ₁ + ... + λ_j) / |λ_{j+1}|`` where ``j`` is the largest
    index with a non-negative cumulative sum (Kaplan & Yorke 1979).

    Parameters
    ----------
    spectrum : array-like
        Lyapunov exponents (any order; sorted descending internally).

    Returns
    -------
    float
        0.0 when every exponent is negative; ``len(spectrum)`` when the
        cumulative sum never turns negative (spectrum incomplete).

    Examples
    --------
    >>> kaplan_yorke_dimension([0.906, 0.0, -14.57])   # Lorenz
    2.062...
    """
    s = np.sort(np.asarray(spectrum, dtype=float))[::-1]
    if s.size == 0 or s[0] < 0.0:
        return 0.0
    cum = np.cumsum(s)
    nonneg = np.nonzero(cum >= 0.0)[0]
    j = int(nonneg[-1])
    if j == s.size - 1:
        return float(s.size)  # spectrum doesn't close — dimension saturates
    return float(j + 1 + cum[j] / abs(s[j + 1]))

Fixed points

fixed_points 

fixed_points(
    map_sys: Any,
    *,
    box: tuple | None = None,
    n_seeds: int = 200,
    tol: float = 1e-12,
    max_iter: int = 60,
    dedup_tol: float = 1e-06,
    seed: int | None = None,
) -> list[FixedPoint]

Find fixed points f(x) = x of a discrete map by multi-start Newton.

Seeds are drawn uniformly from box plus points sampled from a short orbit; each runs Newton on g(x) = f(x) − x with the map's Jacobian. Converged roots are deduplicated and classified by the eigenvalues of J(x*) (stable iff every |λ| < 1).

PARAMETER DESCRIPTION
map_sys

TYPE: DiscreteMap

box

Search box; defaults to the orbit's bounding box padded by 50 %, or [-2, 2]^dim if the orbit diverges.

TYPE: ((lo...), (hi...)) DEFAULT: None

n_seeds

Random seeds (orbit points are added on top).

TYPE: int DEFAULT: 200

tol

Residual tolerance ‖f(x) − x‖.

TYPE: float DEFAULT: 1e-12

max_iter

Newton iterations per seed.

TYPE: int DEFAULT: 60

dedup_tol

Distance below which two roots are considered the same point.

TYPE: float DEFAULT: 1e-06

seed

RNG seed for reproducibility.

TYPE: int DEFAULT: None

RETURNS DESCRIPTION
list[FixedPoint]

Sorted by first coordinate.

Examples:

>>> fps = fixed_points(Henon())
>>> [fp.x[0] for fp in fps]   # analytic: (-0.7 ± sqrt(0.49 + 5.6)) / 2.8
Source code in src/tsdynamics/analysis/fixedpoints/__init__.py 
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
def fixed_points(
    map_sys: Any,
    *,
    box: tuple | None = None,
    n_seeds: int = 200,
    tol: float = 1e-12,
    max_iter: int = 60,
    dedup_tol: float = 1e-6,
    seed: int | None = None,
) -> list[FixedPoint]:
    """
    Find fixed points ``f(x) = x`` of a discrete map by multi-start Newton.

    Seeds are drawn uniformly from ``box`` plus points sampled from a short
    orbit; each runs Newton on ``g(x) = f(x) − x`` with the map's Jacobian.
    Converged roots are deduplicated and classified by the eigenvalues of
    ``J(x*)`` (stable iff every ``|λ| < 1``).

    Parameters
    ----------
    map_sys : DiscreteMap
    box : ((lo...), (hi...)), optional
        Search box; defaults to the orbit's bounding box padded by 50 %,
        or ``[-2, 2]^dim`` if the orbit diverges.
    n_seeds : int
        Random seeds (orbit points are added on top).
    tol : float
        Residual tolerance ``‖f(x) − x‖``.
    max_iter : int
        Newton iterations per seed.
    dedup_tol : float
        Distance below which two roots are considered the same point.
    seed : int, optional
        RNG seed for reproducibility.

    Returns
    -------
    list[FixedPoint]
        Sorted by first coordinate.

    Examples
    --------
    >>> fps = fixed_points(Henon())
    >>> [fp.x[0] for fp in fps]   # analytic: (-0.7 ± sqrt(0.49 + 5.6)) / 2.8
    """
    if not isinstance(map_sys, DiscreteMap):
        raise TypeError("fixed_points currently supports DiscreteMap systems")

    rng = np.random.default_rng(seed)
    dim = map_sys.dim
    params = map_sys.params.as_tuple()
    step_fn = type(map_sys)._step
    jac_fn = type(map_sys)._jacobian

    def f(x: np.ndarray) -> np.ndarray:
        return np.asarray(step_fn(x, *params), dtype=float).ravel()

    def jac_at(x: np.ndarray) -> np.ndarray:
        return np.atleast_2d(np.asarray(jac_fn(x, *params), dtype=float))

    # Seed pool: random box points + on-orbit points
    import contextlib

    orbit_pts = np.empty((0, dim))
    with contextlib.suppress(Exception):  # divergent maps still get box seeds
        orbit_pts = map_sys.copy().iterate(steps=100, max_retries=15).y[20:]

    if box is None:
        if orbit_pts.size:
            lo, hi = orbit_pts.min(axis=0), orbit_pts.max(axis=0)
            span = np.where(hi - lo < 1e-3, 1.0, hi - lo)
            lo, hi = lo - 0.5 * span, hi + 0.5 * span
        else:
            lo, hi = -2.0 * np.ones(dim), 2.0 * np.ones(dim)
    else:
        lo = np.asarray(box[0], dtype=float).reshape(dim)
        hi = np.asarray(box[1], dtype=float).reshape(dim)

    seeds = rng.uniform(lo, hi, size=(n_seeds, dim))
    if orbit_pts.size:
        seeds = np.vstack([seeds, orbit_pts[:: max(1, len(orbit_pts) // 20)]])

    roots: list[np.ndarray] = []
    eye = np.eye(dim)
    for x in seeds:
        x = x.copy()
        ok = False
        for _ in range(max_iter):
            try:
                gx = f(x) - x
            except Exception:  # noqa: BLE001
                break
            if not np.all(np.isfinite(gx)):
                break
            if np.linalg.norm(gx) < tol:
                ok = True
                break
            A = jac_at(x) - eye
            try:
                dx = np.linalg.solve(A, -gx)
            except np.linalg.LinAlgError:
                break
            if not np.all(np.isfinite(dx)) or np.linalg.norm(dx) > 1e6:
                break
            x = x + dx
        if not ok:
            continue
        if not np.all((x >= lo - 1e-6) & (x <= hi + 1e-6)):
            continue  # outside the search box
        if any(np.linalg.norm(x - r) < dedup_tol for r in roots):
            continue
        roots.append(x)

    out = []
    for r in roots:
        eig = np.linalg.eigvals(jac_at(r))
        out.append(FixedPoint(x=r, eigenvalues=eig, stable=bool(np.all(np.abs(eig) < 1.0))))
    out.sort(key=lambda fp: tuple(fp.x))
    return out

FixedPoint dataclass 

FixedPoint(x: ndarray, eigenvalues: ndarray, stable: bool)

A fixed point of a map with its linear stability data.