API¶

The public Quivira interface provides symbolic matrix utilities and helpers for assembling constrained Lagrangian equations of motion. These functions work directly with heyoka.expression objects and NumPy arrays.

quivira.create_block_diagonal(blocks)¶

Create a block-diagonal matrix from square symbolic arrays.

Parameters:: blocks (list) – A sequence of two-dimensional square arrays to place along the main diagonal.
Raises:: ValueError – If a block is not a two-dimensional square array.

Notes

Entries outside the diagonal blocks are filled with heyoka.expression objects representing zero, rather than Python integer zeros. An empty sequence returns an empty (0, 0) object array.

Returns:: An object array containing the block-diagonal matrix.
Return type:: numpy.ndarray

quivira.gaussian_elimination(S, b)¶

Solve a symbolic linear system using sparsity-aware Gaussian elimination.

Solves the linear system

\[S x = b\]

by first computing an LU factorization of the coefficient matrix and then performing forward and backward substitution.

The factorization is computed for each call. To solve multiple systems with the same coefficient matrix more efficiently, use gaussian_factorize() once and reuse the result with gaussian_solve_factored().

Parameters:

S (numpy.ndarray or sequence) – Square coefficient matrix containing symbolic or numeric entries, with shape (n, n).
b (numpy.ndarray or sequence) – Right-hand side vector containing symbolic or numeric entries. Must contain exactly n elements.

Raises:

ValueError – If S is not square or b has an incompatible size.
ZeroDivisionError – If a structurally zero pivot is encountered during factorization.

Notes

The underlying factorization skips known zero entries to reduce unnecessary symbolic operations. No row pivoting is performed, so the solver may fail for invertible matrices with zero pivots.

The input matrix and right-hand side are not modified.

Returns:: Symbolic solution vector with shape (n, 1).
Return type:: numpy.ndarray

quivira.identity(x)¶

Return the input unchanged.

This utility is useful as a no-op placeholder when a function is needed but the caller wants to preserve the original value.

Parameters:: x – Any Python object to return unchanged.
Returns:: The original object passed in.

quivira.invert_matrix(matrix)¶

Return the inverse of a square matrix using Gaussian elimination.

The inverse is defined by

\[A A^{-1} = I,\]

where \(I\) is the identity matrix. The implementation solves the linear systems

\[A x_j = e_j\]

for each basis vector \(e_j\) and stores the resulting columns in the inverse matrix.

Parameters:: matrix (numpy.ndarray or sequence) – Square matrix to invert.
Raises:: ValueError – If matrix is not a square two-dimensional array.

Notes

This routine is intended for the symbolic matrices used throughout Quivira. It relies on the existing Gaussian elimination routine and is therefore compatible with heyoka.expression entries.

Returns:

Inverse matrix with the same object dtype as the: input matrix.

Return type:

numpy.ndarray

Constrained Lagrangian mechanics¶

quivira.build_mass_matrix(lagrangian, qd)¶

Return the symbolic mass matrix for a Lagrangian.

The mass matrix is defined by

\[M_{ij} = \frac{\partial^2 \mathcal L}{\partial \dot q_i \partial \dot q_j}.\]

Parameters:

lagrangian (heyoka.expression) – Scalar Lagrangian \(\mathcal L(q, \dot q)\).
qd (sequence) – Generalized velocities in the same ordering as the model coordinates.

Returns:

Symbolic mass matrix with shape [n, n].

Return type:

numpy.ndarray

quivira.build_Jdqd(J, q, qd)¶

Return the vector \(\dot J \dot q\) from the Jacobian and state rate.

The derivative of the constraint Jacobian along the trajectory is

\[\dot J = \sum_i \frac{\partial J}{\partial q_i} \dot q_i,\]

so the product entering the constrained-acceleration equation is

\[\dot J \dot q = \sum_i \frac{\partial J}{\partial q_i} \dot q_i \dot q.\]

Parameters:

J (numpy.ndarray or sequence) – Constraint Jacobian \(J(q)\).
q (sequence) – Generalized coordinates.
qd (sequence) – Generalized velocities \(\dot q\).

Returns:

The vector \(\dot J \dot q\).

Return type:

numpy.ndarray

quivira.build_constraint_jacobian(F, q)¶

Return the constraint Jacobian \(J = \partial F / \partial q\).

The algebraic constraints are expressed as a vector \(F(q) = 0\). The Jacobian maps coordinate perturbations into constraint-space variations and is used to enforce the acceleration-level constraint \(J \ddot{q} + \dot{J} \dot{q} = 0\).

Parameters:

F (sequence) – Constraint residuals \(F(q)\).
q (sequence) – Generalized coordinates.

Returns:

Jacobian matrix with shape [m, n].

Return type:

numpy.ndarray

quivira.find_constrained_accelerations(qdd_free, J, M, Jdqd, Q_nc=None, *, return_W=False)¶

Return constrained accelerations and the constraint multipliers.

The constrained acceleration is

\[\ddot{q} = \ddot{q}_{\mathrm{free}} - M^{-1} J^T W^{-1} \left(J \ddot{q}_{\mathrm{free}} + \dot{J} \dot{q}\right),\]

where

\[W = J M^{-1} J^T.\]

this is simplified to: The constrained acceleration is

\[\ddot{q} = \ddot{q}_{\mathrm{applied}} + G \lambda,\]

where

\[G = M^{-1} J^T, \qquad W = JG,\]

and

\[W \lambda = -(J \ddot{q}_{\mathrm{applied}} + \dot{J}\dot{q}).\]

The matrix G is calculated by solving M G = J^T, avoiding the explicit inversion of M. The mass matrix is factorized only once, and the factorization is reused for all right-hand sides.

Parameters:

qdd_free (sequence) – Free acceleration vector.
J (numpy.ndarray or sequence) – Constraint Jacobian \(J\).
M (numpy.ndarray or sequence) – Mass matrix \(M\).
Jdqd (sequence) – Product \(\dot J \dot q\).
Q_nc (numpy.ndarray or sequence, optional) – Generalized applied forces, in the same ordering as q. Entries may be symbolic functions of q, qd, and time. Defaults to None.
return_W (bool, optional) – Whether to also return the Delassus matrix \(W\). Default is False.

Returns:

By default, a pair (qdd, lambda_mult) containing the constrained acceleration vector and Lagrange multiplier vector. If return_W is True, a triple (qdd, lambda_mult, W) is returned; \(W = J M^{-1} J^T\) is the Delassus matrix. The multiplier entries follow the constraint-row ordering in \(J\).

Return type:

tuple

quivira.lagrange_eom(lagrangian, F, q, qd, Q_nc=None, *, return_multipliers=False, return_W=False)¶

Assemble the first-order ODE system for the constrained dynamics.

This helper forms the free-acceleration terms from the Lagrangian, evaluates the constraint Jacobian and its time derivative, and then returns the equivalent first-order system that can be passed to a numerical integrator.

Parameters:

lagrangian (heyoka.expression) – Scalar Lagrangian.
F (sequence) – Constraint residuals \(F(q)\).
q (sequence) – Generalized coordinates.
qd (sequence) – Generalized velocities.
Q_nc (sequence, optional) – Generalized applied forces, in the same ordering as q. Entries may be symbolic functions of q, qd, and time. Defaults to None.
return_multipliers (bool, optional) – Whether to also return the Lagrange multiplier vector. Default is False.
return_W (bool, optional) – Whether to also return the Delassus matrix \(W = J M^{-1} J^T\). Default is False.

Returns:

The ODE system as pairs (state_var, rhs_expr). If both options are False, this list is returned directly. If only return_multipliers is True, returns (ode_system, lambda_mult); if only return_W is True, returns (ode_system, W); if both are True, returns (ode_system, lambda_mult, W). The multiplier entries follow the constraint-row ordering in \(J\).

Return type:

list or tuple

Bar joint model class¶

class quivira.bar_joint_model(tolerance=1e-10)¶

Bases: object

Assembly of articulated bars with joints, supports, and external loads.

The model collects the geometric and force data needed to generate a Lagrangian formulation and the constraint equations for the system.

add_bar(length, mass)¶

Append a new bar to the assembly.

Parameters:

length (float) – Bar length.
mass (float) – Bar mass.

add_bar_at_attitude(x0, y0, theta0, length, mass)¶

Append a bar and assign its initial attitude in one step.

Parameters:

x0 (float) – Initial x coordinate of the bar center.
y0 (float) – Initial y coordinate of the bar center.
theta0 (float) – Initial angle of the bar.
length (float) – Bar length.
mass (float) – Bar mass.

add_bar_from_point_pair(start, end, mass_per_length=1.0)¶

Append a bar constructed from its endpoints.

Parameters:

start (tuple) – Cartesian coordinates (x, y) of the first end.
end (tuple) – Cartesian coordinates (x, y) of the second end.
mass_per_length (float, optional) – Distributed mass per unit length. Default is 1.0.

add_external_force(bar_index, s, fx, fy, torque=0.0)¶

Add a point load on a bar.

Parameters:

bar_index (int) – Index of the bar receiving the force.
s (float) – Local coordinate of the force application point.
fx (float) – Force component along the global x axis.
fy (float) – Force component along the global y axis.
torque (float) – Torque mathematically positive

add_joint(bar_index, s)¶

Add a joint linking selected points across bars.

Parameters:

bar_index (list) – Indices of the bars participating in the joint.
s (list) – Local coordinates of the coincident points on each bar.

add_support(bar_index, s, dof_in_x, dof_in_y, x=None, y=None)¶

Add a support constraint on a point of a bar.

Parameters:

bar_index (int) – Index of the constrained bar.
s (float) – Local coordinate of the constrained point.
dof_in_x (bool) – If True, x motion is unconstrained.
dof_in_y (bool) – If True, y motion is unconstrained.
x (float, optional) – Support x position. If omitted, the initial bar position is used. Default is None.
y (float, optional) – Support y position. If omitted, the initial bar position is used. Default is None.

Raises:

ValueError – Raised if the support position is requested before the bar has an initial state.

build_lagrange_eom(return_multipliers=False, return_W=False)¶

Assemble the symbolic equations of motion for the model.

This method builds the system-level Lagrangian, constraint residuals, and external-force terms, then returns the equivalent symbolic first-order ODE system for numerical integration. When return_multipliers or return_W is enabled, the corresponding constraint quantities are included in the return value.

Parameters:

return_multipliers (bool, optional) – If True, also return the symbolic Lagrange multiplier vector associated with the algebraic constraints. Default is False.
return_W (bool, optional) – If True, also return the Delassus matrix \(W = J M^{-1} J^T\), where J is the constraint Jacobian and M is the mass matrix. Default is False.

Returns:

The ODE system as pairs (state_var, rhs_expr). If neither option is enabled, the function returns the ODE system directly. If return_multipliers is True, it returns (ode_system, lambda_mult); if return_W is True, it returns (ode_system, W); and if both are True, it returns (ode_system, lambda_mult, W). The multiplier entries follow the ordering of the constraint residuals in self.F.

Return type:

list or tuple

check_constraint_fulfillment(state)¶

clear()¶: Remove all bars, joints, supports, and forces from the model.

get_initial_state()¶

Return the stacked initial generalized coordinates and velocities.

Returns:: State vector containing [x0, y0, theta0, xd0, yd0, thetad0] for each bar.
Return type:: list

infer_joints()¶

Infer joints from intersections of bars in their initial poses.

Detects intersections anywhere along the bars, groups bars meeting at the same point, and avoids duplicate joints.

Returns:: Total number of joints after inference.
Return type:: int

print_model()¶: Print a human-readable summary of the assembled model.

set_g(g)¶

Set the gravitational acceleration for the model.

Parameters:: g (float) – Gravitational acceleration in the global y direction.

set_init_bar_attitude(bar_index, x0, y0, theta0)¶

Set the initial center position and orientation of a bar.

Parameters:

bar_index (int) – Index of the target bar.
x0 (float) – Initial x coordinate of the bar center.
y0 (float) – Initial y coordinate of the bar center.
theta0 (float) – Initial angle of the bar.

set_init_bar_velocities(bar_index, xd0, yd0, thetad0)¶

Set the initial translational and angular velocities of a bar.

Parameters:

bar_index (int) – Index of the target bar.
xd0 (float) – Initial x velocity of the center.
yd0 (float) – Initial y velocity of the center.
thetad0 (float) – Initial angular velocity.

set_initial_state(initial_state)¶

Set the initial bar poses and velocities from a flattened state vector.

Parameters:: initial_state (list or numpy.ndarray) – State vector in the same order as get_initial_state(), where all generalized coordinates are listed first and all generalized velocities follow.
Raises:: ValueError – Raised when the state length does not match the number of bars or when an individual bar state is malformed.
Returns:: The stacked initial generalized coordinates and velocities.
Return type:: list