Source code for dynamo.vectorfield.FixedPoints

from typing import Optional, Tuple

import numpy as np
from scipy.linalg import eig

[docs]class FixedPoints: """The FixedPoints class stores a list of fixed points and their corresponding Jacobian matrices, and provides methods for computing the eigenvalues of the Jacobian matrices, determining the stability of the fixed points, and identifying saddle/stable fixed points."""
[docs] def __init__(self, X: Optional[np.ndarray] = None, J: Optional[np.ndarray] = None): """This class represents a set of fixed points and their corresponding Jacobian matrices. The fixed points and Jacobian matrices can be provided as arguments, or they can be added to the object later using the add_fixed_points method. The eigvals attribute stores the eigenvalues of the Jacobian matrices, which can be computed using the compute_eigvals method. The is_stable and is_saddle methods can be used to determine the stability and saddle-point status of the fixed points, respectively, and the get_fixed_point_types method returns a list of integers indicating the stability of each fixed point (-1 for stable, 0 for saddle, and 1 for unstable). Args: X: array of fixed points. Defaults to None. J: array of associated jacobians. Defaults to None. """ self.X = X if X is not None else [] self.J = J if J is not None else [] self.eigvals = []
def get_X(self) -> np.ndarray: return np.array(self.X) def get_J(self) -> np.ndarray: return np.array(self.J) def add_fixed_points(self, X: np.ndarray, J: np.ndarray, tol_redundant: float = 1e-4) -> None: for i, x in enumerate(X): redundant = False if tol_redundant is not None and len(self.X) > 0: for y in self.X: if np.linalg.norm(x - y) <= tol_redundant: redundant = True if not redundant: self.X.append(x) self.J.append(J[i]) def compute_eigvals(self) -> None: self.eigvals = [] for i in range(len(self.J)): if self.J[i] is None or np.isnan(self.J[i]).any(): w = np.nan else: w, _ = eig(self.J[i]) self.eigvals.append(w) def is_stable(self) -> np.ndarray: if len(self.eigvals) != len(self.X): self.compute_eigvals() stable = np.ones(len(self.eigvals), dtype=bool) for i, w in enumerate(self.eigvals): if w is None or np.isnan(w).any(): stable[i] = np.nan else: if np.any(np.real(w) >= 0): stable[i] = False return stable def is_saddle(self) -> Tuple[np.ndarray, np.ndarray]: is_stable = self.is_stable() saddle = np.zeros(len(self.eigvals), dtype=bool) for i, w in enumerate(self.eigvals): if w is None or np.isnan(w).any(): saddle[i] = np.nan else: if not is_stable[i] and np.any(np.real(w) < 0): saddle[i] = True return saddle, is_stable def get_fixed_point_types(self) -> np.ndarray: is_saddle, is_stable = self.is_saddle() # -1 -- stable, 0 -- saddle, 1 -- unstable ftype = np.ones(len(self.X)) for i in range(len(ftype)): if self.X[i] is None or np.isnan((self.X[i])).any(): ftype[i] = np.nan else: if is_saddle[i]: ftype[i] = 0 elif is_stable[i]: ftype[i] = -1 return ftype