mathematicskit.ode_dynamics

mathematicskit.ode_dynamics hub · History · Examples · API reference

Fixed-point stability analysis via Jacobian linearization; phase portraits for 2D systems; the logistic map’s period-doubling route to chaos; saddle-node/pitchfork/Hopf bifurcation normal forms; the Van der Pol limit cycle; and Poincare sections of the driven Duffing oscillator.

mathematicskit.ode_dynamics: dynamical systems built on mathematicskit.integrators.

Fixed-point stability analysis via Jacobian linearization and eigenvalue classification (node/saddle/spiral/center); phase portraits for 2D systems; bifurcation diagrams for 1D maps (logistic map, Feigenbaum route to chaos) and canonical parametrized-ODE normal forms (saddle-node, pitchfork, Hopf); limit cycles (Van der Pol oscillator) and Poincare sections for higher-dimensional (periodically driven) flows; Lyapunov’s direct method and Bendixson’s criterion; population (Verhulst logistic, Lotka-Volterra), epidemic (SIR), neuron (FitzHugh-Nagumo), and chemical (Brusselator) models; Kuramoto synchronization; and the Lorenz and Rossler chaotic flows.

class mathematicskit.ode_dynamics.BendixsonResult(X, Y, divergence, rules_out_periodic_orbits=False)

Bases: object

Container for Bendixson’s negative criterion evaluated on a grid.

Parameters:

• X (ndarray)

• Y (ndarray)

• divergence (ndarray)

• rules_out_periodic_orbits (bool)

X: ndarray

Grid x-coordinates.

Type:

ndarray, shape (n, n)

Y: ndarray

Grid y-coordinates.

Type:

ndarray, shape (n, n)

divergence: ndarray

Divergence df/dx + dg/dy of the vector field (f, g).

Type:

ndarray, shape (n, n)

rules_out_periodic_orbits: bool = False

True if the divergence is strictly of one sign on the whole grid, so no closed orbit lies entirely inside the (simply connected) rectangle.

Type:

bool

class mathematicskit.ode_dynamics.Brusselator(state0, a=1.0, b=3.0)

Bases: FlowSystem

The Brusselator reaction scheme.

\(\dot x = a - (b+1)x + x^2y\), \(\dot y = bx - x^2y\), the rate equations of the hypothetical reactions \(A \to X\), \(B + X \to Y + D\), \(2X + Y \to 3X\), \(X \to E\) with the concentrations of \(A, B\) held fixed. The unique fixed point \((a,\ b/a)\) loses stability in a Hopf bifurcation at \(b = 1 + a^2\) (brusselator_hopf_threshold()), beyond which the concentrations oscillate on a stable limit cycle.

Parameters:

• state0 (array-like, shape (2,)) – Initial concentrations (x, y).

• a (float) – Positive feed concentrations.

• b (float) – Positive feed concentrations.

Examples

>>> system = Brusselator([1.0, 1.0], a=1.0, b=1.5)  # b < 1 + a^2 = 2
>>> result = system.integrate((0.0, 100.0), dt=1e-2, method="rk4")
>>> np.round(result.y[-1], 4)
array([1. , 1.5])

fixed_point()

The unique fixed point (a, b/a).

Return type:

ndarray

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.DuffingOscillator(state0, delta=0.3, alpha=-1.0, beta=1.0, gamma=0.3, omega=1.2)

Bases: FlowSystem

Periodically-driven Duffing oscillator.

\(\ddot x + \delta \dot x + \alpha x + \beta x^3 = \gamma \cos(\omega t)\). Explicit time-dependence (through the drive) makes the (x, y) phase portrait alone insufficient to see the full (3D, including drive phase) state space; a stroboscopic Poincare section – sampling (x, y) once per drive period – reduces the long-time behavior back to a 2D picture, revealing periodic points for periodic motion and a fractal attractor for chaotic parameters. See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 12.3.

Parameters:

• state0 (array-like, shape (2,))

• delta (float) – Damping.

• alpha (float) – Linear and cubic stiffness (alpha<0, beta>0 gives the classic double-well potential).

• beta (float) – Linear and cubic stiffness (alpha<0, beta>0 gives the classic double-well potential).

• gamma (float) – Drive amplitude and angular frequency.

• omega (float) – Drive amplitude and angular frequency.

Examples

>>> system = DuffingOscillator([1.0, 0.0], delta=0.3, alpha=-1.0, beta=1.0, gamma=0.5, omega=1.2)
>>> result = system.integrate((0.0, 10.0), dt=1e-3, method="rk4")
>>> result.y.shape[1]
2

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.FitzHughNagumo(state0, a=0.7, b=0.8, eps=0.08, current=0.0)

Bases: FlowSystem

FitzHugh-Nagumo neuron model.

\(\dot v = v - v^3/3 - w + I\), \(\dot w = \varepsilon(v + a - bw)\): a fast “membrane voltage” \(v\) with a cubic nullcline and a slow recovery variable \(w\). For small injected current \(I\) the rest state is stable but excitable – a large enough kick fires a single large spike before returning to rest; for \(I\) between the two values from fitzhugh_nagumo_hopf_currents(), the rest state is unstable and the neuron fires periodically (a relaxation limit cycle).

Parameters:

• state0 (array-like, shape (2,)) – Initial (v, w).

• a (float) – Model parameters (FitzHugh’s classic values 0.7, 0.8, 0.08).

• b (float) – Model parameters (FitzHugh’s classic values 0.7, 0.8, 0.08).

• eps (float) – Model parameters (FitzHugh’s classic values 0.7, 0.8, 0.08).

• current (float) – Injected current \(I\).

Examples

>>> system = FitzHughNagumo([-1.2, -0.6], current=0.0)
>>> result = system.integrate((0.0, 200.0), dt=1e-2, method="rk4")
>>> v_star, w_star = fitzhugh_nagumo_fixed_point(current=0.0)
>>> bool(abs(result.y[-1, 0] - v_star) < 1e-3)
True

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.FixedPointResult(location, jacobian, eigenvalues, classification='', stable=False)

Bases: object

Container for a 2D fixed point’s linear-stability classification.

Parameters:

• location (ndarray)

• jacobian (ndarray)

• eigenvalues (ndarray)

• classification (str)

• stable (bool)

classification: str = ''

One of "stable node", "unstable node", "saddle", "stable spiral", "unstable spiral", "center", "degenerate node".

Type:

str

eigenvalues: ndarray

Eigenvalues of the Jacobian.

Type:

ndarray, shape (2,), complex

jacobian: ndarray

Jacobian evaluated at the fixed point.

Type:

ndarray, shape (2, 2)

location: ndarray

Fixed-point coordinates.

Type:

ndarray, shape (2,)

stable: bool = False

Whether both eigenvalues have negative real part.

Type:

bool

class mathematicskit.ode_dynamics.FlowSystem(state0)

Bases: ABC

Common base for an autonomous flow dy/dt = f(y, t).

Concrete subclasses must set self._rhs_njit (an @njit dispatcher, signature (state, t, params) -> dstate, the mathematicskit.integrators.RHSFunc convention) and self.params in __init__, and call super().__init__(state0).

integrate(t_span, dt=None, method='rk4', **kwargs)

Integrate the flow forward in time.

Parameters:

• t_span (tuple of float) – (t0, t1).

• dt (float, optional) – Fixed step (required for method="rk4"); initial step attempt for method="dopri5" (defaults to (t1-t0)/1000).

• method (str)

• **kwargs – Forwarded to mathematicskit.integrators.dopri5_integrate().

Return type:

OdeTrajectory

Returns:

OdeTrajectory

params: ndarray = array([], dtype=float64)

reset(state0=None, t0=0.0)

Reset state and clock; see e.g. physicskit’s identical pattern.

Parameters:

t0 (float)

abstractmethod rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.KuramotoModel(theta0, omegas, K=1.0)

Bases: FlowSystem

Kuramoto model of \(N\) all-to-all coupled phase oscillators.

\[\dot\theta_i = \omega_i + \frac{K}{N}\sum_{j=1}^N \sin(\theta_j - \theta_i) = \omega_i + K r \sin(\psi - \theta_i),\]

where \(r e^{i\psi} = N^{-1}\sum_j e^{i\theta_j}\) is the complex order parameter (kuramoto_order_parameter()). The mean-field form makes each step \(O(N)\) rather than \(O(N^2)\).

Parameters:

• theta0 (array-like, shape (N,)) – Initial phases.

• omegas (array-like, shape (N,)) – Natural frequencies (passed to the integrator as params).

• K (float) – Coupling strength.

Examples

>>> import numpy as np
>>> system = KuramotoModel([0.0, 1.0, 2.0], omegas=[0.0, 0.0, 0.0], K=1.0)
>>> result = system.integrate((0.0, 50.0), dt=1e-2, method="rk4")
>>> round(float(kuramoto_order_parameter(result.y[-1])[0]), 6)  # identical oscillators lock
1.0

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.Linear2D(state0, A)

Bases: FlowSystem

Linear planar system \(\dot{\mathbf x} = A\mathbf x\).

The textbook object of study for mathematicskit.ode_dynamics.systems.stability.classify_fixed_point_2d(): every trajectory’s qualitative behavior near the origin is fully determined by A’s eigenvalues. See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 5.

Parameters:

• state0 (array-like, shape (2,)) – Initial condition.

• A (array-like, shape (2, 2)) – System matrix.

Examples

>>> import numpy as np
>>> system = Linear2D([1.0, 0.0], A=[[-1.0, 0.0], [0.0, -2.0]])
>>> result = system.integrate((0.0, 5.0), dt=1e-3, method="rk4")
>>> bool(abs(result.y[-1, 0]) < 1e-2)
True

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.LogisticGrowth(state0, r=1.0, K=1.0)

Bases: FlowSystem

Verhulst’s logistic equation \(\dot N = rN(1 - N/K)\).

Growth is exponential while \(N \ll K\) and saturates at the carrying capacity \(K\), a stable fixed point; \(N=0\) is unstable. See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 2.3.

Parameters:

• state0 (array-like, shape (1,)) – Initial population [N0].

• r (float) – Intrinsic growth rate.

• K (float) – Carrying capacity.

Examples

>>> system = LogisticGrowth([10.0], r=1.0, K=100.0)
>>> result = system.integrate((0.0, 20.0), dt=1e-2, method="rk4")
>>> round(float(result.y[-1, 0]), 3)
100.0

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.LogisticMap(r)

Bases: object

The logistic map \(x_{n+1} = r x_n (1 - x_n)\) on [0, 1].

For r in (0, 1) all orbits decay to 0; (1, 3) a single stable fixed point; (3, 3.449...) a stable period-2 cycle, then successive period-doublings accumulating at \(r_\infty \approx 3.56995\), beyond which the map is chaotic (with periodic windows). See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 10.

Parameters:

r (float) – Growth-rate parameter.

Examples

>>> m = LogisticMap(r=2.5)
>>> orbit = m.iterate(x0=0.2, n_transient=200, n_keep=1)
>>> round(float(orbit[0]), 4)
0.6

fixed_points()

Nonzero fixed point \(x^* = 1 - 1/r\) (plus \(x=0\), always a fixed point).

Return type:

ndarray

Returns:

ndarray, shape (2,) – [0.0, 1 - 1/r].

Examples

>>> LogisticMap(r=2.0).fixed_points()
array([0. , 0.5])

iterate(x0, n_transient=500, n_keep=1)

Iterate the map, discarding a transient.

Parameters:

• x0 (float) – Initial condition in (0, 1).

• n_transient (int) – Iterations discarded to let the orbit approach its attractor.

• n_keep (int) – Iterations returned after the transient.

Return type:

ndarray

Returns:

ndarray, shape (n_keep,)

class mathematicskit.ode_dynamics.LorenzSystem(state0, sigma=10.0, rho=28.0, beta=2.6666666666666665)

Bases: FlowSystem

Lorenz’s convection model (Lorenz, 1963).

\(\dot x = \sigma(y - x)\), \(\dot y = x(\rho - z) - y\), \(\dot z = xy - \beta z\). For Lorenz’s classic values \(\sigma = 10\), \(\rho = 28\), \(\beta = 8/3\) every fixed point is unstable and trajectories settle onto the butterfly- shaped strange attractor, where nearby orbits separate exponentially (sensitive dependence on initial conditions).

Parameters:

• state0 (array-like, shape (3,))

• sigma (float)

• rho (float)

• beta (float)

Examples

>>> system = LorenzSystem([1.0, 1.0, 1.0])
>>> result = system.integrate((0.0, 10.0), dt=1e-2, method="rk4")
>>> result.y.shape
(1001, 3)

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.LotkaVolterra(state0, alpha=1.0, beta=1.0, delta=1.0, gamma=1.0)

Bases: FlowSystem

Lotka-Volterra predator-prey system.

\(\dot x = \alpha x - \beta xy\) (prey), \(\dot y = \delta xy - \gamma y\) (predators). The coexistence fixed point \((\gamma/\delta,\ \alpha/\beta)\) is a nonlinear center: every orbit around it is closed, because the quantity returned by lotka_volterra_invariant() is conserved. Small oscillations have period \(2\pi/\sqrt{\alpha\gamma}\). See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 6.4-6.5.

Parameters:

• state0 (array-like, shape (2,)) – Initial (prey, predators).

• alpha (float) – Prey growth, predation, predator conversion, and predator death rates.

• beta (float) – Prey growth, predation, predator conversion, and predator death rates.

• delta (float) – Prey growth, predation, predator conversion, and predator death rates.

• gamma (float) – Prey growth, predation, predator conversion, and predator death rates.

Examples

>>> system = LotkaVolterra([1.0, 0.5], alpha=1.0, beta=1.0, delta=1.0, gamma=1.0)
>>> result = system.integrate((0.0, 10.0), dt=1e-3, method="rk4")
>>> result.y.shape[1]
2
>>> system.fixed_point()
array([1., 1.])

fixed_point()

Coexistence fixed point (gamma/delta, alpha/beta).

Return type:

ndarray

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.LyapunovFunctionResult(P, Q, positive_definite=False)

Bases: object

Container for a quadratic Lyapunov function V(x) = x^T P x.

Parameters:

• P (ndarray)

• Q (ndarray)

• positive_definite (bool)

P: ndarray

Symmetric solution of A^T P + P A = -Q.

Type:

ndarray, shape (n, n)

Q: ndarray

Symmetric positive-definite right-hand side used.

Type:

ndarray, shape (n, n)

positive_definite: bool = False

Whether P is positive definite, i.e. whether V certifies asymptotic stability of the origin (Lyapunov’s theorem).

Type:

bool

class mathematicskit.ode_dynamics.Nonlinear2D(state0, f)

Bases: FlowSystem

General nonlinear planar system \(\dot{\mathbf x} = f(\mathbf x)\).

Wraps an arbitrary plain-Python vector field for integration and phase-portrait visualization (a compiled njit dispatcher is built lazily via numba’s dynamic njit on the supplied callable, so any numba-compatible pure function works).

Parameters:

• state0 (array-like, shape (2,))

• f (Callable[[float, float], tuple]) – f(x, y) -> (dx, dy), numba-compatible (no Python objects).

Examples

>>> import numpy as np
>>> vdp = Nonlinear2D([2.0, 0.0], f=lambda x, y: (y, -x + 0.0 * y))
>>> result = vdp.integrate((0.0, 1.0), dt=1e-3, method="rk4")
>>> result.y.shape[1]
2

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.OdeTrajectory(t, y, method='', extra=<factory>)

Bases: object

Container for the output of a FlowSystem.integrate() call.

Parameters:

• t (ndarray)

• y (ndarray)

• method (str)

• extra (dict)

extra: dict

Free-form diagnostics slot.

Type:

dict

method: str = ''

Integrator used ("rk4" or "dopri5").

Type:

str

t: ndarray

Time samples.

Type:

ndarray, shape (n_steps + 1,)

y: ndarray

State trajectory.

Type:

ndarray, shape (n_steps + 1, dim)

class mathematicskit.ode_dynamics.RosslerSystem(state0, a=0.2, b=0.2, c=5.7)

Bases: FlowSystem

Rossler’s chaotic flow.

\(\dot x = -y - z\), \(\dot y = x + ay\), \(\dot z = b + z(x - c)\). Only one term (\(zx\)) is nonlinear. Trajectories spiral outward in the \((x, y)\) plane until \(x\) exceeds \(c\), when \(z\) spikes and folds them back toward the center – a stretch-and-fold mechanism producing a chaotic attractor for the classic values \(a = b = 0.2\), \(c = 5.7\).

Parameters:

• state0 (array-like, shape (3,))

• a (float)

• b (float)

• c (float)

Examples

>>> system = RosslerSystem([1.0, 1.0, 0.0])
>>> result = system.integrate((0.0, 100.0), dt=1e-2, method="rk4")
>>> result.y.shape
(10001, 3)

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.SIRModel(state0, beta=0.3, gamma=0.1)

Bases: FlowSystem

Kermack-McKendrick SIR model (population fractions).

\(\dot S = -\beta SI\), \(\dot I = \beta SI - \gamma I\), \(\dot R = \gamma I\). The basic reproduction number \(R_0 = \beta/\gamma\) sets the epidemic threshold: infections initially grow only if \(R_0 S_0 > 1\).

Parameters:

• state0 (array-like, shape (3,)) – Initial fractions (S0, I0, R0).

• beta (float) – Transmission rate.

• gamma (float) – Recovery rate.

Examples

>>> system = SIRModel([0.99, 0.01, 0.0], beta=0.3, gamma=0.1)
>>> result = system.integrate((0.0, 100.0), dt=0.1, method="rk4")
>>> round(float(result.y[-1].sum()), 12)  # S + I + R is conserved
1.0

property r0: float

Basic reproduction number beta / gamma.

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

class mathematicskit.ode_dynamics.VanDerPolOscillator(state0, mu=1.0)

Bases: FlowSystem

Van der Pol oscillator \(\ddot x - \mu(1-x^2)\dot x + x = 0\).

Written as a first-order system \(\dot x = y\), \(\dot y = \mu(1-x^2)y - x\). For any \(\mu > 0\), every trajectory except the (unstable) origin converges to a unique stable limit cycle (a consequence of the Poincare-Bendixson theorem, since the system is dissipative for \(|x|>1\) and anti-dissipative for \(|x|<1\)) – unlike a linear center, the cycle’s amplitude and shape are independent of initial conditions. See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 7.4.

Parameters:

• state0 (array-like, shape (2,)) – Initial (x, y).

• mu (float) – Nonlinear damping strength (mu=0 reduces to the harmonic oscillator, a center with no limit cycle).

Examples

>>> system = VanDerPolOscillator([2.0, 0.0], mu=1.0)
>>> result = system.integrate((0.0, 50.0), dt=1e-3, method="rk4")
>>> result.y.shape[1]
2

rhs(state, t=0.0)

Plain-Python evaluation of dy/dt (thin wrapper around self._rhs_njit).

Return type:

ndarray

Parameters:

• state (ndarray)

• t (float)

mathematicskit.ode_dynamics.bendixson_criterion(f, x_range, y_range, n=101, h=1e-06)

Bendixson’s negative criterion on a rectangle.

For a planar system \(\dot x = f(x, y)\), \(\dot y = g(x, y)\), if the divergence \(\partial f/\partial x + \partial g/\partial y\) is not identically zero and does not change sign on a simply connected region \(D\), then no periodic orbit lies entirely in \(D\) (Bendixson, 1901). Proof sketch: by Green’s theorem, the divergence integrated over the region enclosed by a closed orbit equals the flux of the field across the orbit, which is zero because the field is tangent to it. This function samples the divergence (by central differences) on an n x n grid and reports whether it is strictly of one sign there – a numerical check, not a proof.

Parameters:

• f (Callable[[float, float], tuple]) – f(x, y) -> (dx, dy), accepting NumPy arrays elementwise.

• x_range (tuple of float) – (min, max) of the rectangle.

• y_range (tuple of float) – (min, max) of the rectangle.

• n (int) – Grid resolution per axis.

• h (float) – Central-difference step.

Return type:

BendixsonResult

Returns:

BendixsonResult

Examples

>>> damped = lambda x, y: (y, -x - 0.5 * y)  # divergence = -0.5
>>> bendixson_criterion(damped, (-2, 2), (-2, 2)).rules_out_periodic_orbits
True
>>> vdp = lambda x, y: (y, (1 - x**2) * y - x)  # divergence = 1 - x^2
>>> bendixson_criterion(vdp, (-3, 3), (-3, 3)).rules_out_periodic_orbits
False

mathematicskit.ode_dynamics.bifurcation_diagram(r_values, x0=0.5, n_transient=500, n_keep=100)

Compute the logistic map’s bifurcation diagram.

Parameters:

• r_values (array-like of float)

• x0 (float) – Shared initial condition.

• n_transient (int)

• n_keep (int)

Returns:

r_plot, x_plot (ndarray) – Flattened (r, x) pairs suitable for a scatter plot (n_keep points per r value).

Examples

>>> r_plot, x_plot = bifurcation_diagram([2.5], n_transient=200, n_keep=1)
>>> round(float(x_plot[0]), 4)
0.6

mathematicskit.ode_dynamics.brusselator_hopf_threshold(a)

Critical \(b\) for the Brusselator’s Hopf bifurcation.

At \((a, b/a)\) the Jacobian is \(\begin{pmatrix} b - 1 & a^2 \\ -b & -a^2\end{pmatrix}\), with determinant \(a^2 > 0\) and trace \(b - 1 - a^2\). The fixed point is therefore a stable focus/node for \(b < 1 + a^2\) and unstable beyond it, where a limit cycle appears.

Parameters:

a (float)

Return type:

float

Returns:

float – \(b_c = 1 + a^2\).

Examples

>>> brusselator_hopf_threshold(1.0)
2.0

mathematicskit.ode_dynamics.classify_fixed_point_2d(jacobian, location=None)

Classify a 2D linear system’s fixed point from its Jacobian.

For \(\dot{\mathbf x} = J\mathbf x\) (or the linearization of a nonlinear system about a fixed point), let \(\tau = \mathrm{tr}(J)\) and \(\Delta = \det(J)\); the eigenvalues are \(\lambda_{1,2} = (\tau \pm \sqrt{\tau^2 - 4\Delta})/2\). Then:

• \(\Delta < 0\): saddle (always unstable).

• \(\Delta > 0\), \(\tau^2 > 4\Delta\): node (real eigenvalues, same sign as \(\tau\)).

• \(\Delta > 0\), \(\tau^2 < 4\Delta\): spiral (complex eigenvalues, real part sign = sign of \(\tau\)).

• \(\Delta > 0\), \(\tau = 0\): center (purely imaginary).

• \(\tau^2 = 4\Delta \neq 0\): degenerate node (repeated real eigenvalue).

The last two are exact-equality conditions, tested here against a small tolerance rather than against zero. A near-borderline system is genuinely ambiguous rather than merely hard to classify numerically: linearization is inconclusive for centers and degenerate nodes, since an arbitrarily small nonlinear term can push the trajectory either way (Strogatz, Sec. 6.3).

See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 5.2 (the classification diagram in the \((\tau, \Delta)\) plane).

Parameters:

• jacobian (ndarray)

• location (ndarray, shape (2,), optional) – The fixed point’s coordinates, for record-keeping; defaults to the origin.

Return type:

FixedPointResult

Returns:

FixedPointResult

Examples

>>> import numpy as np
>>> result = classify_fixed_point_2d(np.array([[-1.0, 0.0], [0.0, -2.0]]))
>>> result.classification
'stable node'
>>> result2 = classify_fixed_point_2d(np.array([[0.0, 1.0], [-1.0, 0.0]]))
>>> result2.classification
'center'

mathematicskit.ode_dynamics.estimate_feigenbaum_delta(x0=0.5)

Estimate the Feigenbaum constant from the logistic map’s first few period-doubling bifurcation points.

\(\delta = \lim_{n\to\infty} \dfrac{r_n - r_{n-1}}{r_{n+1} - r_n} \approx 4.6692\dots\), universal across a broad class of period-doubling routes to chaos (Feigenbaum, 1978). Using only the first few (numerically located) bifurcation points gives a rough estimate, converging slowly toward the true constant as more bifurcations are included.

Parameters:

x0 (float) – Initial condition used to detect each period’s onset.

Return type:

float

Returns:

float

Examples

>>> delta = estimate_feigenbaum_delta()
>>> 3.0 < delta < 6.0
True

mathematicskit.ode_dynamics.estimate_limit_cycle_amplitude(mu, t_transient=200.0, t_observe=50.0, dt=0.001)

Estimate the Van der Pol limit cycle’s amplitude (max |x|) after transients have decayed.

Parameters:

• mu (float)

• t_transient (float) – Integration time discarded to let the trajectory approach the cycle.

• t_observe (float) – Additional integration time over which the amplitude is measured.

• dt (float)

Return type:

float

Returns:

float – max(|x|) over the observation window.

Examples

>>> amp = estimate_limit_cycle_amplitude(mu=1.0, t_transient=100.0, t_observe=30.0)
>>> 1.8 < amp < 2.2
True

mathematicskit.ode_dynamics.estimate_period(t, x)

Estimate a periodic signal’s period from upward zero crossings.

Finds every t where x crosses zero going from negative to positive (linearly interpolating between samples for sub-step accuracy) and returns the mean spacing between consecutive crossings.

Parameters:

• t (ndarray) – Time samples and the (mean-centered, if necessary) signal.

• x (ndarray) – Time samples and the (mean-centered, if necessary) signal.

Return type:

float

Returns:

float – Estimated period; nan if fewer than 2 crossings are found.

Examples

>>> import numpy as np
>>> t = np.linspace(0.0, 20.0, 20000)
>>> x = np.sin(2.0 * np.pi * t / 3.0)  # period 3.0
>>> round(estimate_period(t, x), 2)
3.0

mathematicskit.ode_dynamics.find_fixed_point_newton(f, x0, tol=1e-10, max_iter=100)

Locate a fixed point of f(x) = 0 via multivariate Newton’s method.

\(x_{k+1} = x_k - J(x_k)^{-1} f(x_k)\), with the linear solve at each step done via mathematicskit.linalg.systems.lu.lu_solve_system() (never an explicit matrix inverse) and the Jacobian obtained from numerical_jacobian(). See Burden & Faires, Numerical Analysis, 10th ed., Ch. 10.2.

Parameters:

• f (Callable[[ndarray], ndarray]) – Vector field f(x) -> ndarray, same shape as x.

• x0 (ndarray) – Initial guess.

• tol (float) – Convergence tolerance on ||x_{k+1} - x_k||.

• max_iter (int)

Returns:

• x (ndarray) – The fixed point found.

• converged (bool)

Examples

>>> import numpy as np
>>> f = lambda x: np.array([x[1] - x[0]**2, x[0] + x[1] - 2.0])
>>> x, converged = find_fixed_point_newton(f, np.array([0.5, 0.5]))
>>> converged
True
>>> np.allclose(f(x), 0.0, atol=1e-8)
True

mathematicskit.ode_dynamics.fitzhugh_nagumo_fixed_point(current=0.0, a=0.7, b=0.8)

The (unique, for \(0 < b < 1\)) fixed point of FitzHugh-Nagumo.

Intersecting the nullclines \(w = v - v^3/3 + I\) and \(w = (v + a)/b\) gives the cubic \(-v^3/3 + (1 - 1/b)v + I - a/b = 0\), which is strictly decreasing (so has one real root) when \(b < 1\).

Parameters:

• current (float)

• a (float)

• b (float)

Return type:

ndarray

Returns:

ndarray, shape (2,) – (v*, w*).

Examples

>>> np.round(fitzhugh_nagumo_fixed_point(0.0), 4)
array([-1.1994, -0.6243])

mathematicskit.ode_dynamics.fitzhugh_nagumo_hopf_currents(a=0.7, b=0.8, eps=0.08)

Injected currents at which the rest state loses/regains stability.

The Jacobian at the fixed point is \(\begin{pmatrix} 1 - v^2 & -1 \\ \varepsilon & -\varepsilon b\end{pmatrix}\), whose determinant \(\varepsilon(1 - b + bv^2)\) is positive for \(b<1\), so stability changes exactly when the trace \(1 - v^2 - \varepsilon b\) vanishes: \(v_H = \pm\sqrt{1 - \varepsilon b}\). Substituting into the nullcline relation gives the Hopf currents

\[I_H = \frac{v_H + a}{b} - v_H + \frac{v_H^3}{3}.\]

Parameters:

• a (float)

• b (float)

• eps (float)

Return type:

ndarray

Returns:

ndarray, shape (2,) – [I_low, I_high]; the rest state is unstable for I_low < I < I_high.

Examples

>>> np.round(fitzhugh_nagumo_hopf_currents(), 4)
array([0.3313, 1.4187])

mathematicskit.ode_dynamics.hopf_limit_cycle_radius(r, kind='supercritical')

Limit-cycle radius of the Hopf normal form (polar-coordinate ODE).

Supercritical Hopf: \(\dot \rho = r\rho - \rho^3\), \(\dot\theta = \omega\); the origin is a stable spiral for \(r<0\) and unstable for \(r>0\), with a stable limit cycle of radius \(\rho = \sqrt r\) born at \(r=0\) and growing as \(\sqrt r\). See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 8.2.

Parameters:

• r (float or array-like of float)

• kind (str) – Only the supercritical case is implemented.

Return type:

ndarray

Returns:

ndarray – sqrt(r) where r > 0, else 0 (no limit cycle).

Examples

>>> hopf_limit_cycle_radius(4.0)
2.0
>>> hopf_limit_cycle_radius(-1.0)
0.0

mathematicskit.ode_dynamics.kuramoto_lorentzian_order_parameter(K, gamma)

Steady-state coherence for Lorentzian natural frequencies (\(N \to \infty\)).

For \(g(\omega) = \gamma / [\pi(\gamma^2 + \omega^2)]\), Kuramoto’s self-consistency equation can be solved exactly: incoherence (\(r = 0\)) is the only state below the critical coupling \(K_c = 2/(\pi g(0)) = 2\gamma\), and above it

\[r = \sqrt{1 - K_c / K}.\]

Parameters:

• K (float or array-like) – Coupling strength.

• gamma (float) – Half-width of the Lorentzian frequency distribution.

Returns:

ndarray or float

Examples

>>> kuramoto_lorentzian_order_parameter(4.0, gamma=1.0)  # K_c = 2
0.7071067811865476
>>> kuramoto_lorentzian_order_parameter(1.0, gamma=1.0)
0.0

mathematicskit.ode_dynamics.kuramoto_order_parameter(theta)

Kuramoto order parameter \(r e^{i\psi} = N^{-1}\sum_j e^{i\theta_j}\).

\(r = 0\) for phases spread uniformly around the circle and \(r = 1\) for perfect phase locking.

Parameters:

theta (array-like, shape (..., N)) – Phases; leading axes (e.g. time) are preserved.

Returns:

r, psi (ndarray, shape (…)) – Coherence and mean phase.

Examples

>>> import numpy as np
>>> r, psi = kuramoto_order_parameter(np.zeros(5))
>>> float(r)
1.0
>>> r, _ = kuramoto_order_parameter(np.linspace(0, 2 * np.pi, 4, endpoint=False))
>>> bool(r < 1e-12)
True

mathematicskit.ode_dynamics.logistic_growth_solution(t, n0, r, K)

Closed-form solution of the logistic equation.

\[N(t) = \frac{K}{1 + \left(\frac{K - N_0}{N_0}\right) e^{-rt}}\]

(Verhulst, 1838), the sigmoid “logistic curve”.

Parameters:

• t (float or array-like) – Time(s).

• n0 (float) – Initial population (> 0).

• r (float) – Growth rate and carrying capacity.

• K (float) – Growth rate and carrying capacity.

Return type:

ndarray

Returns:

ndarray – \(N(t)\).

Examples

>>> float(logistic_growth_solution(0.0, n0=10.0, r=1.0, K=100.0))
10.0
>>> round(float(logistic_growth_solution(50.0, n0=10.0, r=1.0, K=100.0)), 6)
100.0

mathematicskit.ode_dynamics.lorenz_fixed_points(sigma=10.0, rho=28.0, beta=2.6666666666666665)

Closed-form fixed points of the Lorenz system.

The origin is always a fixed point; for \(\rho > 1\) the two convection-roll states

\[C_\pm = \left(\pm\sqrt{\beta(\rho - 1)},\ \pm\sqrt{\beta(\rho - 1)},\ \rho - 1\right)\]

also exist (Lorenz, 1963; Strogatz, Sec. 9.2). sigma does not affect their location, only their stability.

Parameters:

• sigma (float)

• rho (float)

• beta (float)

Return type:

ndarray

Returns:

ndarray, shape (1, 3) or (3, 3) – The origin first, then \(C_+\) and \(C_-\) when \(\rho > 1\).

Examples

>>> fps = lorenz_fixed_points()
>>> fps.shape
(3, 3)
>>> [round(float(v), 4) for v in fps[1]]
[8.4853, 8.4853, 27.0]
>>> lorenz_fixed_points(rho=0.5).shape
(1, 3)

mathematicskit.ode_dynamics.lotka_volterra_invariant(x, y, alpha=1.0, beta=1.0, delta=1.0, gamma=1.0)

Conserved quantity of the Lotka-Volterra system.

\[V(x, y) = \delta x - \gamma \ln x + \beta y - \alpha \ln y\]

Differentiating along a trajectory gives \(\dot V = (\delta - \gamma/x)\dot x + (\beta - \alpha/y)\dot y = 0\), so orbits are level sets of \(V\) – closed curves around the minimum at the coexistence fixed point (Volterra, 1926).

Parameters:

• x (float or array-like) – Prey and predator populations (> 0).

• y (float or array-like) – Prey and predator populations (> 0).

• alpha (float)

• beta (float)

• delta (float)

• gamma (float)

Return type:

ndarray

Returns:

ndarray – \(V(x, y)\).

Examples

>>> float(lotka_volterra_invariant(1.0, 1.0))  # minimum, at the fixed point
2.0

mathematicskit.ode_dynamics.lyapunov_quadratic_form(A, Q=None)

Quadratic Lyapunov function for the linear system \(\dot x = Ax\).

Lyapunov’s direct method (1892) proves stability without solving the ODE: if a function \(V(x) > 0\) (for \(x \neq 0\)) strictly decreases along every trajectory, the origin is asymptotically stable. For a linear system, try \(V(x) = x^T P x\); then \(\dot V = x^T(A^T P + P A)x = -x^T Q x\). Solving the continuous Lyapunov equation

\[A^T P + P A = -Q\]

for a chosen \(Q \succ 0\) gives a positive-definite \(P\) if and only if every eigenvalue of \(A\) has negative real part (Lyapunov’s theorem). The equation is solved with scipy.linalg.solve_continuous_lyapunov() (Bartels-Stewart). See Khalil, Nonlinear Systems, 3rd ed., Thm. 4.6.

Parameters:

• A (array-like, shape (n, n)) – System matrix.

• Q (array-like, shape (n, n), optional) – Symmetric positive-definite matrix; defaults to the identity.

Return type:

LyapunovFunctionResult

Returns:

LyapunovFunctionResult – P, Q, and whether P is positive definite.

Examples

>>> import numpy as np
>>> res = lyapunov_quadratic_form([[-1.0, 0.0], [0.0, -2.0]])
>>> np.round(res.P, 6)
array([[0.5 , 0.  ],
       [0.  , 0.25]])
>>> res.positive_definite
True
>>> lyapunov_quadratic_form([[1.0, 0.0], [0.0, -1.0]]).positive_definite
False

mathematicskit.ode_dynamics.numerical_jacobian(f, x, h=1e-06)

Central-difference Jacobian of a vector field f: R^n -> R^n.

Parameters:

• f (Callable[[ndarray], ndarray])

• x (ndarray)

• h (float) – Step size.

Return type:

ndarray

Returns:

ndarray, shape (n, n)

Examples

>>> import numpy as np
>>> f = lambda x: np.array([x[1], -x[0]])
>>> np.round(numerical_jacobian(f, np.array([0.0, 0.0])), 6)
array([[ 0.,  1.],
       [-1.,  0.]])

mathematicskit.ode_dynamics.pitchfork_fixed_points(r, kind='supercritical')

Fixed points of the pitchfork normal form.

Supercritical: \(\dot x = rx - x^3\), fixed points \(x=0\) (stable for \(r<0\), unstable for \(r>0\)) and, for \(r>0\), \(x=\pm\sqrt r\) (both stable) – one stable branch splitting into two as \(r\) crosses zero. Subcritical: \(\dot x = rx + x^3\), with \(x=\pm\sqrt{-r}\) unstable for \(r<0\). See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 3.4-3.6.

Parameters:

• r (float or array-like of float)

• kind (str)

Return type:

ndarray

Returns:

ndarray, shape (…, 3) – [0, -x*, +x*], with x* = sqrt(|r|) where real, else nan.

Examples

>>> import numpy as np
>>> pitchfork_fixed_points(4.0, kind="supercritical")
array([ 0., -2.,  2.])
>>> bool(np.any(np.isnan(pitchfork_fixed_points(4.0, kind="supercritical"))))
False

mathematicskit.ode_dynamics.rossler_fixed_points(a=0.2, b=0.2, c=5.7)

Closed-form fixed points of the Rossler system.

Setting the right-hand side to zero gives \(z = -y\), \(x = -ay\), and then \(ay^2 + cy + b = 0\), so

\[y_\pm = \frac{-c \pm \sqrt{c^2 - 4ab}}{2a}, \qquad (x, y, z) = (-a y_\pm,\ y_\pm,\ -y_\pm),\]

which are real when \(c^2 \ge 4ab\).

Parameters:

• a (float)

• b (float)

• c (float)

Return type:

ndarray

Returns:

ndarray, shape (2, 3) – One fixed point per row (the one near the origin first).

Examples

>>> fps = rossler_fixed_points()
>>> [round(float(v), 4) for v in fps[0]]
[0.007, -0.0351, 0.0351]
>>> [round(float(v), 4) for v in fps[1]]
[5.693, -28.4649, 28.4649]

mathematicskit.ode_dynamics.saddle_node_fixed_points(r)

Fixed points of the saddle-node normal form \(\dot x = r + x^2\).

Two real fixed points \(x = \pm\sqrt{-r}\) exist for \(r < 0\) (one stable, one unstable), merge at \(r=0\), and disappear for \(r>0\) – the defining “collision and annihilation” signature of a saddle-node bifurcation. See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 3.2.

Parameters:

r (float or array-like of float)

Return type:

ndarray

Returns:

ndarray, shape (…, 2) – [x_stable, x_unstable] (x_stable = -sqrt(-r), x_unstable = +sqrt(-r)); nan where r > 0 (no real fixed points).

Examples

>>> import numpy as np
>>> saddle_node_fixed_points(-4.0)
array([-2.,  2.])
>>> bool(np.all(np.isnan(saddle_node_fixed_points(1.0))))
True

mathematicskit.ode_dynamics.sir_final_size(r0, s0=1.0, i0=0.0)

Fraction of the population still susceptible after the epidemic.

Dividing \(\dot S\) by \(\dot R\) and integrating gives \(S = S_0 e^{-R_0 (R - R_{\text{init}})}\). Since \(I \to 0\), the limit \(S_\infty\) solves the Kermack-McKendrick final-size equation (with \(R_{\text{init}} = 1 - S_0 - I_0\))

\[\ln\frac{S_0}{S_\infty} = R_0\,(S_0 + I_0 - S_\infty),\]

solved here with scipy.optimize.brentq() on \((0, S_0)\).

Parameters:

• r0 (float) – Basic reproduction number beta / gamma.

• s0 (float) – Initial susceptible and infected fractions.

• i0 (float) – Initial susceptible and infected fractions.

Return type:

float

Returns:

float – \(S_\infty\).

Examples

>>> s_inf = sir_final_size(2.0)  # a fully susceptible population
>>> round(s_inf, 4)
0.2032
>>> round(1 - s_inf, 4)  # fraction ever infected
0.7968

mathematicskit.ode_dynamics.sir_peak_infected(r0, s0, i0)

Peak infected fraction of the SIR model.

\(I + S - \ln(S)/R_0\) is conserved, and \(I\) peaks when \(S = 1/R_0\), giving

\[I_{\max} = I_0 + S_0 - \frac{1}{R_0}\left(1 + \ln(R_0 S_0)\right)\]

when \(R_0 S_0 > 1\) (otherwise \(I\) only decreases and \(I_{\max} = I_0\)).

Parameters:

• r0 (float)

• s0 (float) – Initial susceptible and infected fractions.

• i0 (float) – Initial susceptible and infected fractions.

Return type:

float

Returns:

float

Examples

>>> round(sir_peak_infected(3.0, 0.99, 0.01), 4)
0.3038

mathematicskit.ode_dynamics.stroboscopic_poincare_section(system, n_periods, dt=0.001, n_transient_periods=20)

Sample (x, y) once per drive period after discarding transients.

Parameters:

• system (DuffingOscillator)

• n_periods (int) – Number of drive periods to record after the transient.

• dt (float) – Fixed integration step (must divide the drive period reasonably finely).

• n_transient_periods (int) – Number of drive periods discarded before recording begins.

Returns:

xs, ys (ndarray, shape (n_periods,)) – The stroboscopic (once-per-period) samples.

Examples

>>> system = DuffingOscillator([1.0, 0.0], delta=0.3, alpha=-1.0, beta=1.0, gamma=0.3, omega=1.2)
>>> xs, ys = stroboscopic_poincare_section(system, n_periods=5, n_transient_periods=5, dt=1e-2)
>>> len(xs)
5

mathematicskit.ode_dynamics.vector_field_grid(f, x_range, y_range, n=20)

Sample a planar vector field on an n x n grid.

Parameters:

• f (Callable[[float, float], tuple]) – f(x, y) -> (dx, dy) (plain Python; not required to be numba-compatible, unlike Nonlinear2D’s constructor).

• x_range (tuple of float) – (min, max) for each axis.

• y_range (tuple of float) – (min, max) for each axis.

• n (int) – Grid resolution per axis.

Returns:

X, Y, U, V (ndarray, shape (n, n)) – Grid coordinates and vector-field components, suitable for matplotlib.pyplot.quiver(X, Y, U, V).

Examples

>>> X, Y, U, V = vector_field_grid(lambda x, y: (y, -x), (-1, 1), (-1, 1), n=5)
>>> X.shape
(5, 5)