ODEs

Ordinary differential equations are formalized as subclasses of AbstractEquation. The AbstractODE formalizes the syntax, and many concrete subclasses are provided.

class skhippr.odes.AbstractODE.AbstractODE(autonomous: bool, n_dof: int, stability_method=None)

Abstract base class for first-order differential equations. The equilibrium problem can immediately solved bypassing the ODE into the NewtonSolver. If no stability method is provided during instantiation, the default StabilityEquilibrium is used for the equilibrium problem.

Attributes:

autonomousbool

Whether the ODE is autonomous (does not depend on time).

n_dofint

Number of degrees of freedom of the ODE.

xnp.ndarray

State vector.

tfloat

Time variable.

abstractmethod dynamics(t=None, x=None) → ndarray

Return the right-hand side of the first-order ode x_dot = f(t, x) as a 1-D numpy array of size self.n_dof. This method can be passed into scipy.integrate.solve_ivp.

Subclass implementations are expected to check the correct dimensions of x and t using check_dimensions().

Parameters:

• t (float | np.ndarray, optional) – Time variable. If None, self.t is used.

• x (np.ndarray, optional) – State vector. If None, self.x is used.

check_dimensions(t=None, x=None)

Check the dimensions of the time variable and state vector. If t is a scalar, x must be 1-D or have exactly one column. If t is a vector, x must have the same number of columns as t.

Raises:

ValueError – If the dimensions of t and x are incompatible or if x does not have the expected number of rows.

residual_function() → ndarray

Evaluates self.dynamics() at self.t and self.x. The residual is zero at an equilibrium.

stability_criterion(eigenvalues) → bool

Checks if all eigenvalues have a non-positive real part.

derivative(variable, update=False, h_fd=0.0001, t=None, x=None) → ndarray

Provides an interface for optional arguments t and x into the derivative method.

closed_form_derivative(variable, t=None, x=None) → ndarray

Requires subclasses to implement a closed-form derivative with optional arguments t and x.

Nonautonomous ODEs

class skhippr.odes.nonautonomous.Duffing(t: float, x: ndarray, omega: float, alpha: float, beta: float, F: float, delta: float)

Non-autonomous Duffing oscillator as concrete subclass of AbstractODE.

dx[0]/dt = x[1]
dx[1]/dt = -alpha * x[0] - delta * x[1] - beta * [0]**3 + F * cos(omega * t)

Autonomous ODEs

class skhippr.odes.autonomous.Vanderpol(x: ndarray, nu: float, t=0)

Autonomous van der Pol oscillator as a subclass of AbstractODE.

dx[0]/dt = x[1]
dx[1]/dt = nu * (1 - x[0]**2) * x[1] - x[0]

class skhippr.odes.autonomous.Truss(x: ndarray, k: float, c: float, F: float, a: float, l_0: float, m: float)

Autonomous truss system as a subclass of AbstractODE. A mass m can move horizontally with viscous damping (c). It is attached to a linear spring (k, l_0) mounted at the point (0, a). A constant force F acts on the mass. For small values of F, there are three coexisting equilibria. The equations of motion are

dx[0]/dt = x[1]
dx[1]/dt = -k/m * x[0] + k/m * x[0] * l_0 / sqrt(a**2 + x[0]**2) + F/m - c/m * x[1]

class skhippr.odes.autonomous.BlockOnBelt(x: ndarray, epsilon: float, k: float, m: float, Fs: float, vdr: float, delta: float)

Smoothed Block-on-belt system as a subclass of AbstractODE. A block with mass m is placed on a belt with constant velocity vdr. The block is subject to a spring force Fs and smoothed coulomb friction with the belt. The equations of motion are

dx[0]/dt = x[1] dx[1]/dt = -k/m * x[0] + Fs/m * 2/pi * arctan(epsilon * (x[1] - vdr)) / (1 + delta * abs(x[1] - vdr))

Linear time-periodic ODEs

To investigate the properties of the stability methods without influence of a periodic solution, the module skhippr.odes.ltp provides linear time-periodic ordinary differential equations with an equilibrium at zero. Hill stability methods are applicable by searching for the trivial periodic solution X = 0 of the corresponding HBMEquation.

class skhippr.odes.ltp.HillODE(t, x, g_fcn, a=0, b=1, omega=1, damping=0)

Encodes damped Hill-type ODEs of the form

dx[0]/dt = x[1]
dx[1]/dt = - d*x[1] - (a + b*g(omega*t))*x[0]

where g is a 2*pi-periodic function, e.g., a cosine or a rectangular wave.

class skhippr.odes.ltp.HillLTI(t, x, a=0, b=1, damping=0, omega=1)

Subclass of HillODE with g = lambda t: 1, encoding a linear time-invariant (LTI) Hill equation.

fundamental_matrix(t, t_0=None) → ndarray

Compute the fundamental matrix of the linear time-invariant (LTI) Hill system using the matrix exponential.

Parameters:

• t (float) – The final time at which to evaluate the fundamental matrix.

• t_0 (float, optional) – The initial time. If None, defaults to self.t.

class skhippr.odes.ltp.SmoothedMeissner(t, x, smoothing=0, a=0, b=1, omega=1, damping=0)

Subclass of HillODE with g(t) being a smoothed square wave. The smoothing parameter smoothing must lie between 0 and 1.

• For smoothing = 0, we obtain the Meissner equation (Hill equation with rectangular forcing). In this case, a close-form expression for the fundamental matrix is available.

• For smoothing = 1, we obtain the Mathieu equation (Hill equation with cosine forcing)

fundamental_matrix(t_end: float, t_0: float = None) → ndarray

Compute the closed-form fundamental matrix of the Meissner equation (smoothing = 0).

class skhippr.odes.ltp.MathieuODE(t, x, a=0, b=1, omega=1, damping=0)

Subclass of SmoothedMeissner with smoothing = 1. This corresponds to the Mathieu equation, which is a special case of the Hill equation with cosine forcing.

class skhippr.odes.ltp.TruncatedMeissner(t, x, N_harmonics, a=0, b=1, omega=1, damping=0)

Truncated Meissner equation as a subclass of HillODE. The forcing is the Fourier series of the rectangular wave of the Meissner equation, but truncated to the first N_harmonics harmonics.

class skhippr.odes.ltp.ShirleyODE(t, x, E_alpha, E_beta, b, omega)

Two-state quantum system as a subclass of :py:class:`~skhippr.odes.AbstractODE. This is the example discussed by Shirley in https://doi.org/10.1103/PhysRev.138.B979 . The equations of motion are

dx[0]/dt = -i * E_alpha * x[0] - i * 2 * b * cos(omega * t) * x[1]
dx[1]/dt = -i * E_beta * x[1] - i * 2 * b * cos(omega * t) * x[0]

where i is the complex unit.

Note

The equations of motion are complex-valued, so the state vector x is expected to be complex-valued as well. This system can only be handled with the complex-valued HBM formulation.