SpringDamperNode#

Module: maddening.nodes.spring Stability: stable Algorithm ID: MADD-NODE-003 Version: 1.0.0

Summary#

A linear spring with viscous damping connecting a point mass to a moving attachment point. The node owns one end of the spring; the other end arrives each step as the boundary input anchor_position.

Governing Equations#

\[ m \ddot{x} = -k\,(x - a(t) - \ell_0) - c\,\dot{x} \]

where \(x\) is the position of this end of the spring, \(a(t)\) the anchor position, \(\ell_0\) the rest length, \(k\) the stiffness, \(c\) the damping coefficient and \(m\) the point mass.

Discretization#

Semi-implicit (symplectic) Euler: the velocity is advanced first and the position is advanced with the already-updated velocity.

\[ v^{n+1} = v^n + \frac{\Delta t}{m}\left[-k\,(x^n - a^n - \ell_0) - c\,v^n\right], \qquad x^{n+1} = x^n + \Delta t\, v^{n+1} \]

Declared order of accuracy: DiscretizationOrder(spatial=None, temporal=1.0). There is no spatial order — the node integrates an ODE and has no grid. The temporal claim is 1st order globally in both state fields. Unlike a body under a state-independent force, the position does not gain an order here: the spring force depends on the position itself, so the leading error term does not cancel. The claim is measured, not asserted; see Verification Evidence.

Implementation Mapping#

Equation Term	Implementation	Notes
\(-k(x - a - \ell_0)\) (spring force)	`maddening.nodes.spring.SpringDamperNode.update`	`force = -k * (position - anchor - rest)`
\(-c\dot{x}\) (viscous damping)	`maddening.nodes.spring.SpringDamperNode.update`	`- c * velocity`, same expression
\(a(t)\) (anchor position)	`maddening.nodes.spring.SpringDamperNode.boundary_input_spec`	`anchor_position`; defaults to the origin when unsupplied
Velocity update (\(\dot v = F/m\))	`maddening.nodes.spring.SpringDamperNode.update`	`velocity + (force / m) * dt`
Position update (\(\dot x = v\))	`maddening.nodes.spring.SpringDamperNode.update`	`position + velocity * dt`, using the new velocity
Continuous right-hand side	`maddening.nodes.spring.SpringDamperNode.derivatives`	For `maddening.core.simulation.integrators.integrate_node`; not the scheme `update()` uses
Backward-Euler residual	`maddening.nodes.spring.SpringDamperNode.implicit_residual`	\(x^{n+1} - x^n - \Delta t f(x^{n+1})\)
Reaction force delivered over a flux edge	`maddening.nodes.spring.SpringDamperNode.compute_boundary_fluxes`	`spring_force`, same constants as `update`
Coupled-pair limit (warning)	`maddening.nodes.spring._anchored_pair_advisories`	`compile()` warns when two nodes anchored on each other in a coupling group have a converged step that grows, MADD-ANO-098
Converged step of an anchored pair	`maddening.nodes.spring._anchored_pair_step_matrix`	The map on \((x_a - x_b, v_a, v_b)\) the warning takes the spectral radius of

Assumptions and Simplifications#

Linear spring (Hooke’s law) — no hardening or softening
Viscous damping, linear in velocity
Point mass: no rotational dynamics
One translational degree of freedom
No collision detection with other objects

Validated Physical Regimes#

Parameter	Verified Range	Notes
`stiffness`	\(10^{-2}\) – \(10^{6}\) N/m	Very stiff springs need a small \(\Delta t\)
`damping`	\(0\) – \(10^{4}\) N·s/m
`mass`	\(> 0\) kg	Zero mass divides by zero; rejected in the constructor

Known Limitations and Failure Modes#

1st-order integration: the error is \(O(\Delta t)\) in both position and velocity
No stability check: a single node is stable only for \(k\,\Delta t^2 + 2c\,\Delta t < 4m\) (see Stability Conditions), and nothing enforces it at runtime
A calibrated constant reaches derivatives() and implicit_residual() only when it is passed. Both take the injected params by the same {**self.params, **params} rule as update() (MADD-ANO-018, resolved in 0.4.0), so hand the node’s gm.params entry to integrate_node(..., params=) or implicit_euler_step(..., params=); called without it, they integrate the constructor constants
No nonlinear spring behaviour, and no contact
Two nodes anchored on each other in a converged coupling group are unstable for \(k\,\Delta t > c\), far inside the single-node limit, and for \(c\,\Delta t + k\,\Delta t^2 > 2m\). compile() warns about such a pair and does not refuse it (MADD-ANO-098). See Stability Conditions.

Stability Conditions#

One node, fixed or prescribed anchor. The step maps \((x, v)\) by a matrix with trace \(2 - (k\Delta t^2 + c\Delta t)/m\) and determinant \(1 - c\Delta t/m\). The Jury conditions give the exact bound

\[ k\,\Delta t^2 + 2c\,\Delta t < 4m , \]

which is \(\Delta t < 2\sqrt{m/k}\) (\(\omega \Delta t < 2\)) undamped. Damping tightens it.

Two nodes anchored on each other in a converged coupling group (each node’s position is the other’s anchor_position, with rest lengths \(+\ell_0\) and \(-\ell_0\) so that they form one spring). At the converged fixed point each node’s force reads its own position at the old time and its partner’s at the new time. The spring forces therefore do not cancel, and the pair’s momentum is not conserved. For an equal pair (the same \(k\), \(c\) and \(m\)), the sum \(S = v_a + v_b\) decouples exactly:

\[ S^{n+1} = \frac{m - c\,\Delta t}{m - k\,\Delta t^2}\, S^n . \]

The coupling iteration converges for \(k\Delta t^2 < m\); Gauss-Seidel contracts by \((k\Delta t^2/m)^2\) per sweep. Where it converges and \(c\Delta t < m\), the pair is stable only for

\[ c \ge k\,\Delta t , \]

and the factor is exactly 1 at equality. The relative motion \(x_a - x_b\) is stable whenever \(c\,\Delta t < 2m\), so the instability is the pair drifting off together, with a speed that grows every step. The single-node bound does not see it. With \(k = 1000\), \(c = 2\), \(m = 0.5\) and \(\Delta t = 0.01\), the single node has \(k\Delta t^2 + 2c\Delta t = 0.14\) against \(4m = 2\), while the pair’s centre-of-mass velocity grows by a factor of 1.2 per step. Without a coupling group the pair exchanges positions lagged by the schedule. That is a different scheme with its own limit, which this bound does not describe; it also diverges at those constants.

Above \(c\Delta t = m\) the factor is negative, and it passes \(-1\) at \(c\Delta t + k\Delta t^2 = 2m\). So where the iteration converges, an equal pair is stable exactly for

\[ k\,\Delta t \le c \le \frac{2m - k\,\Delta t^2}{\Delta t} . \]

Past the upper figure the pair’s common velocity flips sign and grows every step, though \(c > k\Delta t\): at \(k = 6000\), \(c = 150\), \(m = 1\) and \(\Delta t = 0.01\) each node has \(3.6\) against \(4m = 4\) and \(g = -1.25\). Every equal pair past the single-node bound is past the upper figure too. For an unequal pair the factor has no closed form; the converged step is the linear map on \((x_a - x_b, v_a, v_b)\) that solves the two coupled velocity updates, and its spectral radius is the growth per step. The common translation of the pair, an eigenvalue of exactly 1, is left out of that map.

The compile-time warning. GraphManager.compile() warns, with a UserWarning naming both nodes, their constants, the growth factor and MADD-ANO-098, when it finds two SpringDamperNodes in one coupling group anchored on each other and the spectral radius of their converged step, from the live gm.params values, exceeds 1 by more than one part in a million. It recognises the pattern conservatively: each position reaches the other’s anchor_position through an edge with no transform (or the built-in identity) that is neither additive nor mapped and is the anchor’s only writer; both nodes have one timestep; the group’s max_iterations is at least 2. It does not judge \(k\Delta t^2 \ge m\) (the iteration does not converge, so the converged step is not the one that runs), a wrapped node or a subclass that overrides update, or a value supplied after compile. It warns and never refuses. The full list is in the docstring of maddening.nodes.spring._anchored_pair_advisories.

State Variables#

Field	Shape	Units	Description
`position`	scalar	m	Position of this end of the spring
`velocity`	scalar	m/s	Velocity of this end

Parameters#

Parameter	Type	Default	Units	Description
`stiffness`	float	100.0	N/m	Spring constant \(k\)
`damping`	float	1.0	N·s/m	Damping coefficient \(c\)
`mass`	float	1.0	kg	Point mass \(m\)
`rest_length`	float	1.0	m	Natural length \(\ell_0\)
`initial_position`	float	0.0	m	Starting position
`initial_velocity`	float	0.0	m/s	Starting velocity

Boundary Inputs#

Field	Shape	Default	Description
`anchor_position`	scalar	0.0	Position of the other end of the spring

Boundary Fluxes#

Field	Shape	Units	Description
`spring_force`	scalar	N	Spring-damper force for a downstream node

References#

[@Hairer2006] Hairer, E., Lubich, C. and Wanner, G. (2006). Geometric Numerical Integration. Springer. — Chapter VI derives the order and the energy behaviour of symplectic Euler, which is the scheme this node uses.
[@LeVeque2007] LeVeque, R.J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. — Convergence and stability theory for one-step ODE methods.
[@Roache2002] Roache, P.J. (2002). Code verification by the Method of Manufactured Solutions. J. Fluids Eng. — The method the order claim above is measured with.

Verification Evidence#

Benchmark: MADD-VER-009 — observed temporal order of accuracy by the Method of Manufactured Solutions. A manufactured displacement is injected through anchor_position, which enters the force linearly, so the anchor that makes the trajectory exact is available in closed form. Over a 100/200/400/800 step ladder at fixed final time, in float64, the observed order over the finest pair is 1.029 against the declared 1.0.
Test file: tests/verification/test_mms_order_ode_nodes.py
The same file mutation-tests the study: a mis-scaled timestep, a timestep that drifts with the resolution (a monotone ladder at order 1/2) and a source frozen at \(t=0\) are each required to fail the ladder and name this node.
Stability: tests/nodes/test_spring_stability_limits.py checks both bounds above against the node itself. The single node’s step matrix has spectral radius below 1 just inside \(k\Delta t^2 + 2c\Delta t = 4m\) and above 1 just outside it. The coupled pair’s measured centre-of-mass growth per step equals \((m - c\Delta t)/(m - k\Delta t^2)\) to float32 precision, either side of \(c = k\Delta t\) and at it, and on the heavily damped side. The converged-step matrix the warning judges reproduces 40 steps of the node, for equal and unequal pairs, and the same file checks the warning against the pattern.

Changelog#

Version	Date	Change
1.0.0	2025-03-01	Initial implementation
1.0.0	2026-09-20	Declared order of accuracy added and measured (MADD-VER-009)
1.0.0	2026-10-02	Exact single-node stability bound and the coupled-pair limit \(c \ge k\Delta t\) stated (MADD-ANO-098); numerics unchanged
1.0.0	2026-10-02	The coupled pair’s upper limit \(c\Delta t + k\Delta t^2 \le 2m\) stated, and `compile()` warns about an anchored pair whose converged step grows (MADD-ANO-098); numerics unchanged