physicskit.classical#

physicskit.classical hub · History · Examples · API reference

physicskit.classical: classical_mechanics_kit.

Visual and computational demonstrations of classical mechanics – from Newtonian vector dynamics to Lagrangian variational principles, Hamiltonian phase spaces, coupled lattice chains, and rigid body rotations – built on Numba-accelerated symplectic integrators.

class physicskit.classical.BeadOnRotatingHoop(theta0, thetadot0, R=1.0, omega=2.0, g=9.81)[source]#

Bases: LagrangianSystem

A bead sliding without friction on a hoop of radius R rotating at fixed angular speed Omega about the vertical diameter.

Generalized coordinate: q = (theta,), the bead’s polar angle on the hoop. In the rotating frame the effective Lagrangian is autonomous (no explicit time dependence), so the system is a genuine 1-DOF conservative Hamiltonian; above the critical speed Omega_c = sqrt(g / R) the theta=0 equilibrium becomes unstable and two symmetric stable equilibria appear (a pitchfork bifurcation).

Parameters:

theta0 (float) – Initial angle and angular velocity.
thetadot0 (float) – Initial angle and angular velocity.
R (float) – Hoop radius.
omega (float) – Fixed rotation rate Omega about the vertical diameter.
g (float) – Gravitational acceleration.

energy(state=None)[source]#

Return the total mechanical energy of state.

Parameters:: state (ndarray) – State to evaluate; defaults to self.state.
Return type:: float
Returns:: float – Total mechanical energy.

class physicskit.classical.CoupledOscillators(q0, qdot0, n=3, m=1.0, k=1.0)[source]#

Bases: LagrangianSystem

A chain of N masses connected by linear springs, both ends fixed to immovable walls – the small-N, harmonic prototype of the larger lattice chains in physicskit.classical.systems.chains.

Generalized coordinates: q_i = displacement of mass i from equilibrium, i = 1..N.

Parameters:

q0 (array-like, shape (n,)) – Initial displacements and velocities.
qdot0 (array-like, shape (n,)) – Initial displacements and velocities.
n (int) – Number of masses.
m (float) – Mass and spring constant (uniform across the chain).
k (float) – Mass and spring constant (uniform across the chain).

energy(state=None)[source]#

Return the total mechanical energy of state.

Parameters:: state (ndarray) – State to evaluate; defaults to self.state.
Return type:: float
Returns:: float – Total mechanical energy.

class physicskit.classical.DoublePendulum(theta0, thetadot0, m1=1.0, m2=1.0, l1=1.0, l2=1.0, g=9.81)[source]#

Bases: LagrangianSystem

Planar double pendulum: two point masses on massless rods.

Generalized coordinates: q = (theta1, theta2), angles from the downward vertical. A classic showcase of deterministic chaos and of the non-separable-Hamiltonian case for symplectic integration (the mass matrix depends on theta2 - theta1).

Parameters:

theta0 (array-like, shape (2,)) – Initial angles and angular velocities.
thetadot0 (array-like, shape (2,)) – Initial angles and angular velocities.
m1 (float) – Bob masses.
m2 (float) – Bob masses.
l1 (float) – Rod lengths.
l2 (float) – Rod lengths.
g (float) – Gravitational acceleration.

energy(state=None)[source]#

Return the total mechanical energy of state.

Parameters:: state (ndarray) – State to evaluate; defaults to self.state.
Return type:: float
Returns:: float – Total mechanical energy.

positions(q=None)[source]#

Cartesian (x1, y1, x2, y2) bob positions for a given (theta1, theta2).

Parameters:: q (ndarray) – (theta1, theta2); defaults to the current state.
Returns:: ndarray, shape (4,) – (x1, y1, x2, y2).

class physicskit.classical.DynamicalSystem(state0)[source]#

Bases: ABC

Common base for any system with a state vector evolving in time.

abstractmethod derivatives(t, state)[source]#

Return dstate/dt at (t, state).

Parameters:

t (float) – Current time.
state (ndarray) – Current state vector.

Return type:

ndarray

Returns:

ndarray – Time derivative of state.

abstractmethod energy(state=None)[source]#

Return the total mechanical energy of state.

Parameters:: state (Optional[ndarray]) – State to evaluate; defaults to self.state.
Return type:: float
Returns:: float – Total mechanical energy.

abstractmethod integrate(t_span, dt, method='rk4')[source]#

Integrate the system forward in time.

Parameters:

t_span (tuple of float) – (t0, t1), start and end time.
dt (float) – Fixed step size.
method (str) – Integrator to use; valid values depend on the subclass.

Return type:

SimulationResult

Returns:

SimulationResult – The full trajectory and diagnostics.

reset(state0=None, t0=0.0)[source]#

Reset the system’s state and clock.

Parameters:

state0 (array-like, optional) – New state; if omitted, the current state is kept.
t0 (float) – New time.

Returns:

ndarray – The (possibly updated) current state.

class physicskit.classical.EulerTop(omega0, I1=1.0, I2=2.0, I3=3.0, quat0=None)[source]#

Bases: ODESystem

Torque-free rigid body: Euler’s equations for omega = (w1, w2, w3) in the body frame, coupled to the orientation quaternion (qw, qx, qy, qz).

State: [w1, w2, w3, qw, qx, qy, qz]. The default integrator, implicit_midpoint, exactly conserves both the rotational kinetic energy and the (squared) angular momentum magnitude and quaternion norm – all quadratic invariants of this quadratic ODE (a classical result: the implicit midpoint / average-vector-field map preserves every quadratic invariant of any ODE, not only linear ones).

Excite a spin dominantly about the intermediate-inertia axis (e.g. omega0 = [0.01, 1.0, 0.01] with I1 < I2 < I3) to see the periodic tumbling of the Intermediate Axis Theorem.

Parameters:

omega0 (array-like, shape (3,)) – Initial body-frame angular velocity.
I1 (float) – Principal moments of inertia.
I2 (float) – Principal moments of inertia.
I3 (float) – Principal moments of inertia.
quat0 (array-like, shape (4,), optional) – Initial orientation quaternion (w, x, y, z); defaults to the identity orientation. Normalized automatically.

angular_momentum_squared(state=None)[source]#

|L|^2 in the body frame.

Parameters:: state (ndarray) – State to evaluate; defaults to the current state.
Return type:: float
Returns:: float

energy(state=None)[source]#

Rotational kinetic energy.

Parameters:: state (ndarray) – State to evaluate; defaults to the current state.
Return type:: float
Returns:: float

property omega: ndarray#

property quaternion: ndarray#

rotation_matrix(state=None)[source]#

3x3 body-to-world rotation matrix from the current quaternion.

Parameters:: state (ndarray) – State to evaluate; defaults to the current state.
Return type:: ndarray
Returns:: ndarray, shape (3, 3)

class physicskit.classical.FPUTChain(n=32, m=1.0, k=1.0, beta=0.0, mode=1, amplitude=1.0)[source]#

Bases: _ChainBase

Fermi-Pasta-Ulam-Tsingou beta-lattice: N masses (default 32) with the quartic non-linear coupling potential V(r) = (k/2) r^2 + (beta/4) r^4, r = q_{i+1} - q_i.

Famous for its “paradox”: exciting a single low-order normal mode does not thermalize energy across all modes as naive ergodic reasoning predicts; instead the energy in modal_energies() returns almost exactly to the initially excited mode after a characteristic recurrence time (the FPUT recurrence).

Parameters:

n (int) – Number of masses.
m (float) – Mass and (linear) spring constant.
k (float) – Mass and (linear) spring constant.
beta (float) – Quartic non-linearity coefficient.
mode (int) – Which normal mode to excite initially.
amplitude (float) – Initial amplitude of that mode.

potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

class physicskit.classical.HamiltonianSystem(q0, p0)[source]#

Bases: DynamicalSystem

Base class for canonical Hamiltonian systems H(q, p, t).

Subclasses set ndof, separable, and mass_inv, and build self._force_njit(q, t) -> array (required when separable is True, giving -dV/dq for H = T(p) + V(q)) and/or self._deriv_njit(t, y) -> dy for the full stacked state (required for non-separable systems, and optional otherwise – if omitted it is built automatically from _force_njit).

Parameters:

q0 (array-like, shape (ndof,)) – Initial position and momentum.
p0 (array-like, shape (ndof,)) – Initial position and momentum.

Notes

state is stored flattened as [q_1..q_n, p_1..p_n].

derivatives(t, state)[source]#

Return dstate/dt at (t, state).

Parameters:

t (float) – Current time.
state (ndarray) – Current state vector.

Return type:

ndarray

Returns:

ndarray – Time derivative of state.

energy(state=None)[source]#

Return the total mechanical energy of state.

Parameters:: state (Optional[ndarray]) – State to evaluate; defaults to self.state.
Return type:: float
Returns:: float – Total mechanical energy.

force(q, t=0.0)[source]#

Return -dV/dq (only defined for separable systems).

Parameters:

q (ndarray) – Generalized positions.
t (float) – Current time.

Return type:

ndarray

Returns:

ndarray

Raises:

NotImplementedError – If the system is non-separable (no _force_njit set).

integrate(t_span, dt, method='yoshida4')[source]#

Integrate the system forward in time.

Parameters:

t_span (tuple of float) – (t0, t1), start and end time.
dt (float) – Fixed step size.
method (str) – Integrator to use. "yoshida4"/"verlet" require self.separable.

Return type:

SimulationResult

Returns:

SimulationResult – The full (q, p) trajectory and energy diagnostics.

abstractmethod kinetic_energy(p)[source]#

Return T(p).

Parameters:: p (ndarray) – Canonical momenta.
Return type:: float
Returns:: float

mass_inv = 1.0#

1/m per coordinate; dq/dt = mass_inv * p.

Type:: float or ndarray

ndof: int = 0#

Number of degrees of freedom, inferred from q0.

Type:: int

property p: ndarray#

Canonical momenta (the last ndof entries of state).

Type:: ndarray

abstractmethod potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

property q: ndarray#

Generalized positions (the first ndof entries of state).

Type:: ndarray

separable: bool = True#

Whether H = T(p) + V(q), enabling the Verlet/Yoshida4 integrators.

Type:: bool

static split(state, ndof)[source]#

Split a stacked [q, p] state vector into its two halves.

Parameters:

state (ndarray)
ndof (int)

Returns:

q, p (ndarray, shape (ndof,))

class physicskit.classical.HarmonicChain(q0, p0, m=1.0, k=1.0)[source]#

Bases: _ChainBase

Linear chain of N masses coupled by identical springs k, both ends fixed to walls – the exactly-solvable reference lattice whose normal modes are the sine waves used to analyze the non-linear FPUTChain below.

Parameters:

q0 (array-like, shape (n,)) – Initial displacements and momenta.
p0 (array-like, shape (n,)) – Initial displacements and momenta.
m (float) – Mass and spring constant (uniform across the chain).
k (float) – Mass and spring constant (uniform across the chain).

potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

class physicskit.classical.HeavySymmetricTop(angles0, angledots0, I1=1.0, I3=0.5, M=1.0, l=1.0, g=9.81)[source]#

Bases: LagrangianSystem

Symmetric top with one point fixed, spinning under gravity.

Generalized coordinates q = (phi, theta, psi) – the standard z-x-z Euler angles (precession, nutation, spin). Lagrangian:

L = (I1/2)(thetadot^2 + phidot^2 sin^2(theta))

(I3/2)(psidot + phidot cos(theta))^2 - M*g*l*cos(theta)

where I1 is the transverse moment of inertia about the fixed point, I3 the axial moment, M the mass, and l the distance from the pivot to the center of mass. Because M(q) depends on theta, this is a genuinely non-separable Hamiltonian system – exactly the case physicskit.classical.utils.symbolic.LagrangianEngine and implicit_midpoint were built to handle.

Parameters:

angles0 (array-like, shape (3,)) – Initial (phi, theta, psi) and their time derivatives.
angledots0 (array-like, shape (3,)) – Initial (phi, theta, psi) and their time derivatives.
I1 (float) – Transverse moment of inertia about the fixed point.
I3 (float) – Axial moment of inertia.
M (float) – Mass.
l (float) – Distance from the pivot to the center of mass.
g (float) – Gravitational acceleration.

energy(state=None)[source]#

Return the total mechanical energy of state.

Parameters:: state (ndarray) – State to evaluate; defaults to self.state.
Return type:: float
Returns:: float – Total mechanical energy.

precession_nutation_rate(state=None)[source]#

Instantaneous (phidot, thetadot) – precession and nutation rates.

Parameters:: state (ndarray) – State to evaluate; defaults to the current state.
Returns:: phidot, thetadot (float)

class physicskit.classical.HenonHeilesSystem(q0, p0)[source]#

Bases: HamiltonianSystem

H = (px^2 + py^2)/2 + (x^2 + y^2)/2 + x^2*y - y^3/3.

A non-integrable 2-DOF oscillator (originally a model of stellar orbits in an axisymmetric galactic potential). Below E ~ 1/6 the Poincare section (x, px) at y=0 is dominated by smooth KAM tori; as E approaches and exceeds 1/6 the tori progressively break up into chaotic seas.

Parameters:

q0 (array-like, shape (2,)) – Initial position (x, y) and momentum (px, py).
p0 (array-like, shape (2,)) – Initial position (x, y) and momentum (px, py).

kinetic_energy(p)[source]#

Return T(p).

Parameters:: p (ndarray) – Canonical momenta.
Return type:: float
Returns:: float

mass_inv = 1.0#

1/m per coordinate; dq/dt = mass_inv * p.

Type:: float or ndarray

potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

separable: bool = True#

Whether H = T(p) + V(q), enabling the Verlet/Yoshida4 integrators.

Type:: bool

class physicskit.classical.KeplerSystem(q0, p0, k=1.0, mu=1.0, eps=0.0, c_eps=0.0, c_pn=0.0, softening=0.0)[source]#

Bases: HamiltonianSystem

Planar 2-body Kepler problem with optional non-Newtonian perturbations.

Parameters:

q0 (array-like, shape (2,)) – Initial position and (reduced-mass) momentum.
p0 (array-like, shape (2,)) – Initial position and (reduced-mass) momentum.
k (float) – Gravitational coupling G * M * mu (standard Kepler strength).
mu (float) – Reduced mass (sets kinetic energy T = |p|^2 / (2 mu)).
eps (float) – Exponent and coefficient of an additional c_eps / r**(1+eps) perturbing potential (eps=0 reduces to an extra inverse-square term).
c_eps (float) – Exponent and coefficient of an additional c_eps / r**(1+eps) perturbing potential (eps=0 reduces to an extra inverse-square term).
c_pn (float) – Coefficient of the effective post-Newtonian c_pn / r**3 correction driving perihelion precession.
softening (float) – Small Plummer-style softening length to avoid a singular force at r=0 (0 for the exact point-mass limit).

angular_momentum(q=None, p=None)[source]#

Out-of-plane angular momentum L = x*py - y*px.

Parameters:

q (ndarray) – State to evaluate; defaults to the current state.
p (ndarray) – State to evaluate; defaults to the current state.

Return type:

float

Returns:

float

classmethod from_orbital_elements(a, e, k=1.0, mu=1.0, **kwargs)[source]#

Construct a system starting at perihelion (pure two-body reference orbit).

Parameters:

a (float) – Semi-major axis.
e (float) – Eccentricity.
k (float) – Gravitational coupling and reduced mass, as in __init__.
mu (float) – Gravitational coupling and reduced mass, as in __init__.
**kwargs – Forwarded to __init__ (e.g. eps, c_eps, c_pn, softening).

Returns:

KeplerSystem

kinetic_energy(p)[source]#

Return T(p).

Parameters:: p (ndarray) – Canonical momenta.
Return type:: float
Returns:: float

lrl_vector(q=None, p=None)[source]#

Laplace-Runge-Lenz vector A = p x L - mu*k*r_hat (planar form).

Exactly conserved (a fixed vector pointing at perihelion) for the unperturbed 1/r potential; its slow rotation once c_eps or c_pn is nonzero is the precession signature (e.g. Mercury’s perihelion advance).

Parameters:

q (ndarray) – State to evaluate; defaults to the current state.
p (ndarray) – State to evaluate; defaults to the current state.

Return type:

ndarray

Returns:

ndarray, shape (2,)

potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

separable: bool = True#

Whether H = T(p) + V(q), enabling the Verlet/Yoshida4 integrators.

Type:: bool

class physicskit.classical.LagrangianSystem(q0, qdot0)[source]#

Bases: DynamicalSystem

Base class for systems defined via a Lagrangian L(q, qdot, t).

Subclasses typically obtain their equations of motion from physicskit.classical.utils.symbolic.LagrangianEngine, which derives both a direct acceleration function (for reporting/plotting in (q, qdot) coordinates, via self._accel_njit) and a canonical Hamiltonian form (self._canonical_deriv_njit, used for symplectic integration via implicit_midpoint – valid even when the mass matrix M(q) is configuration-dependent and the system is therefore non-separable, e.g. the double pendulum).

Parameters:

q0 (array-like, shape (ndof,)) – Initial generalized position and velocity.
qdot0 (array-like, shape (ndof,)) – Initial generalized position and velocity.

Notes

state is stored flattened as [q_1..q_n, qdot_1..qdot_n].

acceleration(q, qdot, t)[source]#

Return qddot at (q, qdot, t) via the Euler-Lagrange equations.

Parameters:

q (ndarray)
qdot (ndarray)
t (float)

Return type:

ndarray

Returns:

ndarray

derivatives(t, state)[source]#

Return dstate/dt at (t, state).

Parameters:

t (float) – Current time.
state (ndarray) – Current state vector.

Return type:

ndarray

Returns:

ndarray – Time derivative of state.

abstractmethod energy(state=None)[source]#

Return the total mechanical energy of state.

Parameters:: state (Optional[ndarray]) – State to evaluate; defaults to self.state.
Return type:: float
Returns:: float – Total mechanical energy.

integrate(t_span, dt, method='implicit_midpoint')[source]#

Integrate the system forward in time.

Parameters:

t_span (tuple of float) – (t0, t1), start and end time.
dt (float) – Fixed step size.
method (str) – Integrator to use. "implicit_midpoint" (the default) is symplectic even for non-separable systems, via the Legendre-transformed canonical form; "rk4" integrates the direct (q, qdot) acceleration form and is not symplectic.

Return type:

SimulationResult

Returns:

SimulationResult – The full (q, qdot) trajectory and energy diagnostics (with result.p holding qdot, not a canonical momentum).

momentum(q, qdot)[source]#

Return the canonical momentum p = dL/dqdot at (q, qdot).

Parameters:

q (ndarray)
qdot (ndarray)

Return type:

ndarray

Returns:

ndarray

ndof: int = 0#

Number of degrees of freedom, inferred from q0.

Type:: int

property q: ndarray#

Generalized positions (the first ndof entries of state).

Type:: ndarray

property qdot: ndarray#

Generalized velocities (the last ndof entries of state).

Type:: ndarray

static split(state, ndof)[source]#

Split a stacked [q, qdot] state vector into its two halves.

Parameters:

state (ndarray)
ndof (int)

Returns:

q, qdot (ndarray, shape (ndof,))

class physicskit.classical.ODESystem(state0)[source]#

Bases: DynamicalSystem

A general (possibly non-conservative) first-order ODE system.

Concrete subclasses must set self._deriv_njit (an @njit dispatcher with signature (t, state) -> dstate) in __init__. Suitable for systems that are not naturally Hamiltonian/symplectic, e.g. central-force problems with drag, or rigid-body Euler equations (whose implicit-midpoint flow exactly conserves the quadratic energy and angular-momentum-squared invariants).

derivatives(t, state)[source]#

Return dstate/dt at (t, state).

Parameters:

t (float) – Current time.
state (ndarray) – Current state vector.

Return type:

ndarray

Returns:

ndarray – Time derivative of state.

integrate(t_span, dt, method='implicit_midpoint')[source]#

Integrate the system forward in time.

Parameters:

t_span (tuple of float) – (t0, t1), start and end time.
dt (float) – Fixed step size.
method (str) – Integrator to use.

Return type:

SimulationResult

Returns:

SimulationResult – The full trajectory and energy diagnostics.

class physicskit.classical.PendulumSwarm(q0, p0, g_over_l=1.0)[source]#

Bases: HamiltonianSystem

An ensemble of N independent, identical simple pendulums, H_i = p_i^2 / 2 - (g/l) cos(q_i) for each i = 1..N.

The pendulums are mutually decoupled (each is its own 1-DOF Hamiltonian system) but are integrated together as one N-DOF separable Hamiltonian so a whole ensemble can be evolved with a single Yoshida4/Verlet call. Seeding N points inside a small (dq, dp) box and watching the occupied phase-space patch shear over time – its area invariant even as its shape stretches into a filament – is the standard visual proof of Liouville’s theorem.

Parameters:

q0 (array-like, shape (n,)) – Initial angles and momenta, one pair per pendulum in the ensemble.
p0 (array-like, shape (n,)) – Initial angles and momenta, one pair per pendulum in the ensemble.
g_over_l (float) – Ratio g/l shared by every pendulum.

classmethod from_box(q_center, p_center, dq, dp, n=1000, g_over_l=1.0, seed=0)[source]#

Seed N phase points uniformly inside a (dq x dp) box centered at (q_center, p_center).

Parameters:

q_center (float) – Box center.
p_center (float) – Box center.
dq (float) – Box width and height.
dp (float) – Box width and height.
n (int) – Number of phase points to seed.
g_over_l (float) – Ratio g/l shared by every pendulum.
seed (int) – Random seed for reproducibility.

Returns:

PendulumSwarm

kinetic_energy(p)[source]#

Return T(p).

Parameters:: p (ndarray) – Canonical momenta.
Return type:: float
Returns:: float

per_particle_energy(q=None, p=None)[source]#

Energy of each individual pendulum, used to color/track the swarm.

Conserved per-particle since the ensemble is decoupled.

Parameters:

q (ndarray) – State to evaluate; defaults to the current state.
p (ndarray) – State to evaluate; defaults to the current state.

Return type:

ndarray

Returns:

ndarray, shape (n,)

phase_space_area(q=None, p=None)[source]#

Convex-hull area of the current swarm’s (q, p) point cloud.

A practical proxy for the occupied phase-space volume.

Parameters:

q (ndarray) – State to evaluate; defaults to the current state.
p (ndarray) – State to evaluate; defaults to the current state.

Return type:

float

Returns:

float

potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

separable: bool = True#

Whether H = T(p) + V(q), enabling the Verlet/Yoshida4 integrators.

Type:: bool

class physicskit.classical.ProjectileMotion(q0, p0, m=1.0, g=9.81, drag_coeff=0.0)[source]#

Bases: HamiltonianSystem

Free fall under constant gravity with an oblique initial velocity – the “cannonball” problem – with optional quadratic air drag.

Without drag this is the separable Hamiltonian H(q, p) = |p|^2/(2m) + m*g*y, q = (x, y), whose exact solution is Galileo’s parabolic trajectory (see analytic_trajectory()); because the force is constant, the symplectic Verlet/Yoshida4 backends reproduce it to machine precision at any step size.

With drag_coeff > 0 a quadratic air-resistance force F = -drag_coeff * |v| * v is added. This makes the system genuinely dissipative (energy decreases monotonically, there is no closed-form trajectory), so it is only available through the rk4/implicit_midpoint methods – verlet/yoshida4 would silently ignore it (their force callback only sees the conservative gravity term) and are therefore refused outright.

Parameters:

q0 (array-like, shape (2,)) – Initial position and momentum, q = (x, y), p = (px, py).
p0 (array-like, shape (2,)) – Initial position and momentum, q = (x, y), p = (px, py).
m (float) – Projectile mass and gravitational acceleration.
g (float) – Projectile mass and gravitational acceleration.
drag_coeff (float) – Quadratic air-drag coefficient (0 for the vacuum/exact case).

static analytic_trajectory(speed, angle_deg, g, t, height=0.0)[source]#

Exact vacuum trajectory x(t), y(t) – Galileo’s parabola.

Used to validate the numerical integration (drag_coeff=0 only).

Parameters:

speed (float) – Launch speed and angle above the horizontal, in degrees.
angle_deg (float) – Launch speed and angle above the horizontal, in degrees.
g (float) – Gravitational acceleration.
t (array-like) – Times at which to evaluate the trajectory.
height (float) – Initial height.

Returns:

x, y (ndarray)

classmethod from_launch(speed, angle_deg, m=1.0, g=9.81, drag_coeff=0.0, height=0.0)[source]#

Construct a system launched from (0, height) – the natural “aim the cannon” entry point.

Parameters:

speed (float) – Launch speed.
angle_deg (float) – Launch angle above the horizontal, in degrees.
m (float) – Projectile mass and gravitational acceleration.
g (float) – Projectile mass and gravitational acceleration.
drag_coeff (float) – Quadratic air-drag coefficient (0 for the vacuum/exact case).
height (float) – Initial height.

Returns:

ProjectileMotion

integrate(t_span, dt, method='yoshida4')[source]#

Integrate the system forward in time.

Parameters:

t_span (tuple of float) – (t0, t1), start and end time.
dt (float) – Fixed step size.
method (str) – Integrator to use. "yoshida4"/"verlet" require self.separable.

Returns:

SimulationResult – The full (q, p) trajectory and energy diagnostics.

kinetic_energy(p)[source]#

Return T(p).

Parameters:: p (ndarray) – Canonical momenta.
Return type:: float
Returns:: float

potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

static range_and_max_height(speed, angle_deg, g, height=0.0)[source]#

Closed-form (vacuum) horizontal range and maximum height.

Parameters:

speed (float) – Launch speed and angle above the horizontal, in degrees.
angle_deg (float) – Launch speed and angle above the horizontal, in degrees.
g (float) – Gravitational acceleration.
height (float) – Initial height.

Returns:

range, max_height (float)

separable: bool = True#

Whether H = T(p) + V(q), enabling the Verlet/Yoshida4 integrators.

Type:: bool

class physicskit.classical.SimulationResult(t, y, q=None, p=None, energy=None, method='', extra=<factory>)[source]#

Bases: object

Container for the output of a system’s integrate call.

Parameters:

t (ndarray)
y (ndarray)
q (ndarray | None)
p (ndarray | None)
energy (ndarray | None)
method (str)
extra (dict)

energy: ndarray | None = None#

Total mechanical energy H(t) (or T+V) at each sample.

Type:: ndarray, shape (n_steps + 1,), optional

extra: dict#

Free-form slot for any additional diagnostics a system chooses to attach.

Type:: dict

method: str = ''#

Name of the integrator used (e.g. "yoshida4").

Type:: str

p: ndarray | None = None#

Canonical momenta for Hamiltonian systems, or qdot for Lagrangian systems.

Type:: ndarray, shape (n_steps + 1, ndof), optional

q: ndarray | None = None#

Generalized positions, for Hamiltonian/Lagrangian systems.

Type:: ndarray, shape (n_steps + 1, ndof), optional

t: ndarray#

Time samples.

Type:: ndarray, shape (n_steps + 1,)

y: ndarray#

The raw stacked-state trajectory, in whatever layout the producing system uses internally.

Type:: ndarray, shape (n_steps + 1, state_dim)

class physicskit.classical.SineGordonChain(n=200, m=1.0, k=1.0, q0=None, p0=None)[source]#

Bases: _ChainBase

Discretized pendulum chain (Frenkel-Kontorova model): N coupled pendulums, V = sum_i [k/2 (q_{i+1}-q_i)^2 + (1 - cos q_i)], with free (Neumann-like) boundary conditions so a topological kink can propagate off either end without an artificial restoring wall.

Continuum-limit kink initial data from kink() (q_i = 4*atan(exp((i - i0)/width)), width=1 for the default m=1, k=1) launches a soliton that propagates and scatters off other kinks/antikinks – the discrete analogue of the exact sine-Gordon soliton solutions.

Parameters:

n (int) – Number of pendulums.
m (float) – Mass and coupling constant.
k (float) – Mass and coupling constant.
q0 (array-like, shape (n,), optional) – Initial displacements and momenta; defaults to all zeros. Use kink() to build soliton initial data instead.
p0 (array-like, shape (n,), optional) – Initial displacements and momenta; defaults to all zeros. Use kink() to build soliton initial data instead.

static kink(n, center, width=1.0, velocity=0.0, polarity=1)[source]#

Continuum-limit single-kink (or antikink, polarity=-1) initial condition, optionally boosted to move at velocity (sites per unit time; to the right for positive velocity, assuming the default m=1, k=1 – i.e. a natural wave/”light” speed c = sqrt(k/m) = 1, the sine-Gordon equation’s exact Lorentz-invariance scale, so |velocity| < 1 is required).

width is not a free shape knob: q_xx = sin(q) (the continuum equation this chain approximates, with the default m=1, k=1) is solved by q=4*atan(exp(x)) only for exactly unit width. Passing a different width here without also scaling k to match (the static kink of m q_tt = k q_xx - sin(q) has width = sqrt(k); m only sets the wave speed c = sqrt(k/m)) gives a profile that visibly relaxes/radiates under the true dynamics instead of propagating as a clean, stable soliton.

The exact traveling-wave solution of the continuum sine-Gordon equation is the Lorentz-contracted profile q(x, t) = polarity * 4*atan(exp(gamma * (x - center - v*t) / width)), gamma = 1 / sqrt(1 - v**2) – both the spatial profile at t=0 and its time derivative must be boosted together (not just the momentum on top of the unboosted static profile, which is only a small-v approximation and becomes badly wrong as |v| approaches 1).

Parameters:

n (int) – Chain length.
center (float) – Lattice site the kink is centered on at t=0.
width (float) – Soliton width; must equal sqrt(k) (1, for the default k=1) for a genuine soliton solution.
velocity (float) – Boost velocity, |velocity| < 1.
polarity (int) – +1 for a kink, -1 for an antikink.

Return type:

tuple

Returns:

q0, p0 (ndarray, shape (n,)) – Initial displacements and momenta. p0 is dq/dt at t=0, i.e. the momentum for the default m=1; multiply by m for another mass.

Raises:

ValueError – If abs(velocity) >= 1.

potential_energy(q)[source]#

Return V(q).

Parameters:: q (ndarray) – Generalized positions.
Return type:: float
Returns:: float

physicskit.classical.effective_potential_symmetric_top(theta, p_phi, p_psi, I1, I3, Mgl)[source]#

V_eff(theta) for the heavy symmetric top, given the two conserved (cyclic-coordinate) momenta p_phi and p_psi:

V_eff(theta) = (p_phi - p_psi cos(theta))^2 / (2 I1 sin^2(theta))

p_psi^2 / (2 I3) + Mgl * cos(theta)

Motion in theta is confined to where the total energy E >= V_eff(theta); the turning points bound the nutation range.

Parameters:

theta (array-like) – Angle(s) to evaluate at.
p_phi (float) – Conserved (cyclic-coordinate) momenta.
p_psi (float) – Conserved (cyclic-coordinate) momenta.
I1 (float) – Transverse and axial moments of inertia.
I3 (float) – Transverse and axial moments of inertia.
Mgl (float) – M * g * l (mass times gravity times pivot-to-CM distance).

Returns:

ndarray

physicskit.classical.find_theta_equilibrium(p_phi, p_psi, I1, I3, Mgl, bracket=(0.001, 3.1405926535897932))[source]#

The nutation angle theta_eq at which V_eff(theta) is minimized.

The angle steady (non-nutating) precession would hold at, and the center about which a real top’s nutation oscillates.

Parameters:

p_phi (float) – Conserved (cyclic-coordinate) momenta.
p_psi (float) – Conserved (cyclic-coordinate) momenta.
I1 (float) – Transverse and axial moments of inertia.
I3 (float) – Transverse and axial moments of inertia.
Mgl (float) – M * g * l.
bracket (tuple of float) – Search interval for the minimization.

Return type:

float

Returns:

float

physicskit.classical.nutation_frequency(p_phi, p_psi, I1, I3, Mgl, theta_eq=None, dtheta=1e-05)[source]#

Small-oscillation nutation frequency about theta_eq.

Treating theta as a 1-DOF particle of “mass” I1 moving in the effective potential V_eff (see effective_potential_symmetric_top()), the harmonic-approximation frequency about its minimum is sqrt(V_eff''(theta_eq) / I1), found here by a central finite difference. If theta_eq is not supplied it is located with find_theta_equilibrium().

Parameters:

p_phi (float) – Conserved (cyclic-coordinate) momenta.
p_psi (float) – Conserved (cyclic-coordinate) momenta.
I1 (float) – Transverse and axial moments of inertia.
I3 (float) – Transverse and axial moments of inertia.
Mgl (float) – M * g * l.
theta_eq (float) – Equilibrium angle; computed via find_theta_equilibrium() if omitted.
dtheta (float) – Finite-difference step.

Return type:

float

Returns:

float

physicskit.classical.pendulum_action_angle(E, g_over_l=1.0)[source]#

Action J(E) and period T(E) for a simple pendulum H = p^2/2 - (g/l)cos(q).

For librating orbits (-g/l < E < g/l), uses the standard elliptic-integral closed form: with k^2 = (E + g/l) / (2 g/l),

J(E) = (8/pi) * sqrt(g/l) * [E_ellip(k) - (1 - k^2) * K_ellip(k)] T(E) = 4 * K_ellip(k) / sqrt(g/l)

where K, E are the complete elliptic integrals of the first/second kind (SciPy’s special.ellipk/ellipe take the parameter m = k^2, not the modulus k).

Parameters:

E (float) – Energy, with -g_over_l < E < g_over_l (librating orbits only).
g_over_l (float) – Ratio g/l.

Returns:

J (float) – Action.
T (float) – Orbital period.

Raises:

ValueError – If E is outside the librating range.

physicskit.classical.precession_frequency(p_phi, p_psi, I1, theta_eq)[source]#

Mean precession rate phidot at nutation angle theta_eq: (p_phi - p_psi*cos(theta_eq)) / (I1*sin(theta_eq)**2) – the standard steady-precession-rate approximation, exact when the top is not nutating at all (theta held fixed at theta_eq).

Parameters:

p_phi (float) – Conserved (cyclic-coordinate) momenta.
p_psi (float) – Conserved (cyclic-coordinate) momenta.
I1 (float) – Transverse moment of inertia.
theta_eq (float) – Nutation angle to evaluate at.

Return type:

float

Returns:

float