Zaki Harari
Novelant Scientific Research Inc., Toronto, Canada.
DOI: 10.4236/ijmnta.2026.153007   PDF    HTML   XML   4 Downloads   60 Views   Citations

Abstract

We present the Geometric Quadrature Method (GQM), a coordinate-free, spatial-domain framework for the dynamics of a particle constrained to an arbitrary smooth planar curve. By adopting arc-length parametrization and the Frenet–Serret frame, GQM eliminates the coordinate singularities inherent in Cartesian systems, which fail at vertical tangents or multi-valued trajectory regions. The method operates as a three-pass pipeline. Pass 1 delivers the spatial speed field and normal contact force in exact closed form for any curve with an analytic height function—no quadrature, no approximation, including the full rail loading map and liftoff conditions. Pass 2 confines the time-of-flight computation, which is irreducibly non-elementary by Liouville’s theorem, to a single adaptive quadrature evaluated to machine precision. Pass 3 recovers phase-stable long-time trajectories for conservative systems by exact half-period tiling, bounding accumulated phase error many orders of magnitude below standard explicit or symplectic ODE integrators. The pipeline is anchored on a speed-field kernel encoding the force model in modular form: redefining the kernel alone accommodates friction, quadratic drag, and rotating-frame potentials, leaving the quadrature pipeline intact. Six such kernels are derived. Numerical validation includes a simple pendulum benchmarked against RK4, with symplectic (Störmer-Verlet) behaviour discussed analytically via its known error scaling, and a cubic curve with an inflection point and Coulomb friction, demonstrating that the continuously signed Frenet frame handles curvature sign changes without frame-flip artefacts; the latter result is cross-checked against an independent direct time-domain integration of the same friction model. A unified table reduces eight standard curve families to this single framework.

Keywords

Constrained Dynamics, Arc-Length Parametrization, Frenet-Serret Frame, Nonlinear Dynamics, Geometric Quadrature, Speed-Field Kernel, Inflection Points, Coulomb Friction Kernel

Share and Cite:

1. Introduction

The motion of a particle constrained to a smooth curve under gravity has been a cornerstone of classical mechanics since the seventeenth century, motivating the brachistochrone problem [1], the calculus of variations [2], and Lagrangian mechanics [3]. The brachistochrone itself has spurred a rich line of research into optimal descent curves under dissipation: Šalinić [4] obtained a closed-form solution for descent with Coulomb friction using variational calculus, and Šalinić et al. [5] subsequently extended this to arbitrary conservative and non-conservative force fields via optimal control. More recently, Cherkasov et al. [6] unified the brachistochrone and the Goddard thrust-maximisation problem for a mass point subject to gravity and viscous friction, constructing the optimal synthesis in slope-angle-velocity-mass space. These variational treatments all require the normal contact force as an ingredient of the friction term, making the efficient computation of $N\left(s\right)$ as a function of arc length a prerequisite for, and not merely a by-product of, the trajectory problem.

Despite this maturity, computational treatments for arbitrary planar curves remain fragmented. A persistent bottleneck is the use of the horizontal coordinate $x$ as the independent variable, which introduces the geometric factor $\sqrt{1+{{\phi }^{\prime }}^2}$ into the equations of motion. This term diverges at vertical tangents, fails for closed curves such as ellipses, and obscures the underlying differential geometry—a deficiency long recognised in multibody dynamics [7] but not systematically resolved for the single-particle constrained-motion setting.

The Frenet-Serret frame offers a natural remedy. Shabana [7] studied curvature singularities arising in Frenet-frame formulations and introduced Frenet-oscillation concepts for general motion-trajectory analysis. Bettamin et al. [8] applied the Frenet frame to decompose contact forces along recorded railroad vehicle trajectories, demonstrating the practical value of frame-based force analysis. Friction in such curved-contact settings is governed by the normal force; Marques et al. [9] provide a comprehensive comparison of Coulomb and viscous friction models for multibody systems, establishing the theoretical foundation that the GQM’s kernel table (Section 3) draws upon. For the numerical treatment of constrained Hamiltonian systems over long times, Hairer et al. [10] establish that symplectic integrators conserve a modified (shadow) Hamiltonian, suppressing secular energy growth but leaving period errors of $O\left(h^p\right)$ per cycle, so that phase drift still accumulates as $O\left(h^p{T}_{\text{final}}\right)$ —precisely the limitation that the GQM’s exact half-period tiling is designed to overcome. In astrodynamics, the Binet equation [11] and Kepler-equation inversion [12] provide classical examples of spatial-domain techniques that bypass temporal integration entirely, and the GQM’s Pass-2 time-of-flight quadrature is the direct analogue for general planar curves.

The present work unifies these perspectives into a robust, coordinate-invariant framework built on three interlocking contributions:

• Three-Pass Pipeline: Separation into 1) exact spatial quantities via the speed-field identity, 2) a single scalar quadrature for time-of-flight, and 3) exact half-period tiling for conservative (periodic) systems yields phase stability unattainable by direct ODE integration. The tiling strategy, analogous to Kepler-equation inversion [12], bounds accumulated phase error at $N{\epsilon }_{\tau }$ regardless of the number of periods, whereas symplectic methods reduce energy drift but do not eliminate phase drift [10].

• Speed-Field Identity: The squared speed $V^2\left(s\right)=V_0^2-2g\left[z\left(s\right)-z_0\right]$ is identified as the primary physical descriptor. Combined with the Frenet-Serret normal projection, it yields the full rail loading map $N\left(s\right)$ in exact closed form—a result not available from prior arc-length treatments [7] [8], and a prerequisite for the friction kernels of [4] [9].

• Unified Kernel Table: Speed-field kernels are derived for six physical regimes: frictionless sliding, rigid rolling, Coulomb friction, quadratic drag, surface-tension adhesion, and non-inertial rotating frames. Each kernel is a modular component in the three-pass pipeline, making extension to new force laws straightforward. Together with an eight-family dimensional reduction table, this provides a single reusable framework.

The paper is organized as follows. Table 1 establishes the notation. Section 2 develops the arc-length formulation and the three-pass pipeline, preceded by Section 2.1 which traces the method’s conceptual origin. Section 3 derives the speed-field kernels. Section 4 presents the dimensional reduction table. Section 5 provides numerical validation, followed by a discussion in Section 6 and conclusions in Section 7.

Table 1. Table of notation.

Symbol	Definition	Unit
$s$	Arc-length parameter along the curve	m
$\lambda$	Generic curve parameter (angle, polar angle, etc.)	varies
$r\left(s\right)$	Position vector $\left\{x\left(s\right),z\left(s\right)\right\}$	m
$x\left(s\right),z\left(s\right)$	Horizontal and vertical coordinates	m
$\widehat{t}\left(s\right)$	Unit tangent vector $\text{d}r/\text{d}s$	-
$\widehat{n}\left(s\right)$	Unit normal vector $\left(-\text{d}z/\text{d}s,\text{d}x/\text{d}s\right)$	-
$\kappa \left(s\right)$	Signed curvature	m−1
$\alpha \left(s\right)$	Inclination angle $=\text{arctan}\left(\text{d}z/\text{d}x\right)$	rad
$V\left(s\right)$	Arc-length speed $\text{d}s/\text{d}t$	m∙s−1
$V_0$	Initial speed at $s=s_0$	m∙s−1
$z_0$	Initial height $z\left(s_0\right)$	m
$s_{\star }$	Turning-point arc-length where $V=0$	m
$N\left(s\right)$	Normal contact force (positive = inward)	N
$m$	Particle mass	kg
$g$	Gravitational acceleration	m∙s−2
$\mu$	Coulomb friction coefficient	-
$\beta$	Rolling inertia ratio $I/\left(m{R}^2\right)$	-
$b$	Quadratic drag coefficient ($F_{\text{drag}}=-b{V}^2$ )	kg∙m−1
$\gamma$	Surface tension coefficient (droplet model)	N m−1
${\ell }_c$	Characteristic contact length (droplet model)	m
$\mathcal{K}\left(s\right)$	Speed-field kernel (influence kernel)	m∙s−2
$\mathcal{A}\left(s\right)$	Generalized speed-field kernel	m∙s−2
${\mathcal{K}}_g\left(s\right)$	Gravitational speed-field kernel, $=-g\text{d}z/\text{d}s$	m∙s−2
$T$	Oscillation period	s
$\text{d}s/\text{d}\lambda$	Arc-length element	m

2. Mathematical Formulation

2.1. Conceptual Origin of the Method

We developed the GQM as a spatial-domain framework by building on the simple case of constrained motion: a mass sliding on a flat inclined plane with a constant slope $\alpha$ . In this configuration, Newton’s second law projected onto the surface tangent yields a constant acceleration $a_s=-g\sin \alpha$ , resulting in an exact closed-form quadratic trajectory. This solution requires no numerical approximation, as it follows entirely from the geometry of the slope and the initial conditions. We show that the GQM successfully preserves this exactness even when the slope varies, extending the precision of the inclined plane to arbitrary, non-linear curves.

2.1.1. Stitching Planes: Discrete Convolution in Space

For a piecewise-linear surface with $N$ segments of slope ${\alpha }_k$ , the solution applies exactly on each segment, and the exit state of segment $k$ becomes the entry state of segment $k+1$ through a kinematic handoff. The full trajectory is a superposition of $N$ parabolic arcs gated by Heaviside functions:

$x\left(t\right)={\displaystyle \sum }_{k=0}^{N-1}\left[x_k^{*}+{\dot{x}}_k^{*}\left(t-t_k^{*}\right)-\frac{g}{2}\sin {\alpha }_k\cos {\alpha }_k{\left(t-t_k^{*}\right)}^2\right]\left[H\left(t-t_k^{*}\right)-H\left(t-t_{k+1}^{*}\right)\right].$ (1)

This is a discrete spatial convolution: each slope element contributes an acceleration impulse that propagates forward over the residual time $\left(t-t_k^{\text{*}}\right)$ , weighted by the causal gate. The convolution kernel $\left(t-\tau \right)$ is the Green’s function of ${\text{d}}^2/\text{d}{t}^2$ , applied here in the spatial domain.

2.1.2. Continuous Limit and the Kernel Structure

As $\text{Δ}{x}_k\to \text{d}x$ , (1) passes to a spatial integral. Using the tangential equation of motion $V\text{d}V/\text{d}s=-g\text{d}z/\text{d}s$ and multiplying both sides by $\text{d}s/\text{d}x=\sqrt{1+{{\phi }^{\prime }}^2}$ gives $V\text{d}V/\text{d}x=-g{\phi }^{\prime }\left(x\right)$ , which is a first-order ODE in $x$ for the arc-length speed $V=\text{d}s/\text{d}t$ . Integrating from $x_0$ to $x$ yields the spatial speed field and arrival time:

$V\left(x\right)=\sqrt{V_0^2-2g{\displaystyle {\int }_{x_0}^x{\phi }^{\prime }\left(u\right)\text{d}u}}=\sqrt{V_0^2-2g\left[\phi \left(x\right)-z_0\right]},t\left(x\right)={\displaystyle {\int }_{x_0}^x\frac{\text{d}\xi }{V\left(\xi \right)}}.$ (2)

The arc-length framework in Section 2.5 generalizes (2) to any smooth plane curve via $V^2\left(s\right)=V_0^2-2g\left[z\left(s\right)-z_0\right]$ .

The integrand of (2) encapsulates the method’s design principle. The speed-field kernel

$\mathcal{K}\left(x\right)=g{\phi }^{\prime }\left(x\right)=g\tan \alpha \left(x\right)$ (3)

factors the force law ($g$ ) from the geometry (${\phi }^{\prime }$ ) multiplicatively. The kernel $\mathcal{K}$ acts as an influence kernel: it encodes how a gravitational increment at location $x$ contributes to the cumulative speed change downstream, in exact analogy with the discrete convolution kernel of (1). Replacing gravity with a different physical force—Coulomb friction, quadratic drag, a rotating-frame potential—modifies only $\mathcal{K}$ , leaving the quadrature pipeline intact. The kernel is thus a plug-and-play module, a property formalized and tabulated in Section 3. Equation (2) is the work-energy theorem expressed in Cartesian form: the kernel $g{\phi }^{\prime }\left(x\right)=g\text{d}z/\text{d}x$ is the rate of height gain, and the integral $g{\displaystyle \int {\phi }^{\prime }\text{d}x}=g\left[\phi \left(x\right)-z_0\right]$ is precisely the gravitational potential increment. The arc-length formulation of Section 2.5 generalizes this directly: $-g\text{d}z/\text{d}s$ in (12) is the same kernel expressed in arc-length coordinates.

2.2. Arc-Length Parametrization and the Frenet-Serret Frame

Let $\mathcal{C}$ be a smooth planar curve represented by a unit-speed parametrization

$r:\left[s_0,s_1\right]\to {ℝ}^2,r\left(s\right)=\left(x\left(s\right),z\left(s\right)\right),$ (4)

where $s$ is arc length, so that the unit-speed identity

${\left(\frac{\text{d}x}{\text{d}s}\right)}^2+{\left(\frac{\text{d}z}{\text{d}s}\right)}^2=1$ (5)

holds by definition. This single identity replaces all occurrences of $\sqrt{1+{{\phi }^{\prime }}^2}$ in Cartesian formulations and is valid regardless of whether $\mathcal{C}$ is a graph, a closed curve, or has vertical tangents. The Frenet-Serret frame consists of

$\widehat{t}\left(s\right)=\frac{\text{d}r}{\text{d}s}=\left(\frac{\text{d}x}{\text{d}s},\frac{\text{d}z}{\text{d}s}\right),$(6)

$\widehat{n}\left(s\right)=\left(-\frac{\text{d}z}{\text{d}s},\frac{\text{d}x}{\text{d}s}\right),$(7)

with $\widehat{t}\cdot \widehat{n}=0$ and $\left|\widehat{t}\right|=\left|\widehat{n}\right|=1$ by (5). The signed curvature is defined as

$\kappa \left(s\right)=\frac{\text{d}x}{\text{d}s}\frac{{\text{d}}^{2}z}{\text{d}{s}^2}-\frac{\text{d}z}{\text{d}s}\frac{{\text{d}}^{2}x}{\text{d}{s}^2},$ (8)

and the Frenet-Serret equations read $\text{d}\widehat{t}/\text{d}s=\kappa \widehat{n}$ and $\text{d}\widehat{n}/\text{d}s=-\kappa \widehat{t}$.

Remark (Globally oriented normal). In the classical Frenet frame, $\widehat{n}$ is required to point toward the center of curvature, forcing a 180˚ discontinuous flip whenever the curve passes through an inflection point (where $\kappa =0$ ). GQM avoids this by defining $\widehat{n}$ in (7) as a fixed counter-clockwise rotation of the tangent, independent of the sign of $\kappa$ . This globally oriented normal remains continuous and well-defined through inflection points. The signed curvature $\kappa \left(s\right)$ then carries the full geometric information: positive $\kappa$ means the curve bends in the $\widehat{n}$ direction, negative $\kappa$ means it bends opposite. As a result, the normal contact force $N\left(s\right)$ in (15) transitions continuously between inward and outward loading without any special treatment at the inflection boundary (see Example 2 and Section 3.2).

Remark. When $\mathcal{C}$ is a Cartesian graph $z=\phi \left(x\right)$ , one has $\text{d}x/\text{d}s={\left(1+{{\phi }^{\prime }}^2\right)}^{-1/2}$ , $\text{d}z/\text{d}s={\phi }^{\prime }{\left(1+{{\phi }^{\prime }}^2\right)}^{-1/2}$ , and $\kappa ={\phi }^{″}{\left(1+{{\phi }^{\prime }}^2\right)}^{-3/2}$ , recovering the standard Cartesian expressions as a special case.

2.3. Scope and Standing Assumptions

Before proceeding, we state explicitly the conditions under which each pass of the pipeline (Section 2.7) is valid. These assumptions are used silently throughout Sections 3 - 5 and are collected here for reference.

(A1) Regularity and smoothness of the curve.

The curve $\mathcal{C}$ is assumed at least $C^2$ on the domain of interest, i.e. $x\left(\lambda \right)$ and $z\left(\lambda \right)$ possess continuous second derivatives, so that the signed curvature (8) is well-defined and continuous. Pass-1 closed-form results (Theorem 1, Theorem 3) require only $z\left(\lambda \right)$ to be differentiable; the $C^2$ requirement is needed specifically for $\kappa \left(\lambda \right)$ and hence for the normal force (15). Curves with $C^1$ but not $C^2$ regularity (e.g. a cusp or a corner where the tangent direction is discontinuous) fall outside the present formulation; see the extended discussion in Section 6.7.

(A2) Existence of an arc-length parametrization.

The map $\lambda \mapsto r\left(\lambda \right)$ is assumed regular, i.e.$\left|\text{d}r/\text{d}\lambda \right|\ne 0$ on the domain of interest, so that arc length $s\left(\lambda \right)={\displaystyle \int \left|\text{d}r/\text{d}{\lambda }^{\prime }\right|\text{d}{\lambda }^{\prime }}$ is a strictly increasing, and hence invertible, function of $\lambda$ . This guarantees that the unit-speed parametrization (4) exists and that $s$ may be used as the independent variable throughout Section 2.7. The curve is further assumed rectifiable (finite arc length on any bounded sub-interval), which follows automatically from $C^1$ regularity on a compact domain.

(A3) Conservative vs. dissipative force models.

Theorem 1, Proposition 2, and the exact half-period tiling of Pass 3 (18) all assume that the only tangential force is the conservative gravitational projection $-g\text{d}z/\text{d}s$ (or, more generally, the gradient of a fixed potential $U\left(x,z\right)$ , as in the rotating-frame kernel of Section 3). Under this assumption the motion is time-reversible and periodic between turning points, which is what permits tiling. When a dissipative term is present (Coulomb friction, quadratic drag, or the droplet kernel of Table 2), the speed-squared field $V^2\left(s\right)$ is no longer a single-valued function of position alone—it depends on the direction and history of motion through $\text{sign}\left(V\right)$ and $\text{sign}\left(N\right)$ —so the motion is in general non-periodic and Pass 3 reduces to direct interpolation on the Pass-2 table without tiling, exactly as already noted after (18) and in Section 6.6.

Table 2. Speed-field kernels $\mathcal{A}\left(s\right)$ for six physical models. $\widehat{t}$ and $\widehat{n}$ are unit tangent and normal; $\kappa$ is signed curvature; $\Theta \left(s\right)={\displaystyle {\int }_{s_0}^s\kappa \text{d}{s}^{\prime }}$ is the total turning angle. The viscous drag row uses quadratic (Rayleigh) drag $F=-b{V}^2$ ; ${\mathcal{K}}_g\left(s\right)=-g\text{d}z/\text{d}s$ .

Model	Kernel $\mathcal{A}\left(s\right)$	$V^2\left(s\right)$
Frictionless sliding	$-g\frac{\text{d}z}{\text{d}s}$	$V_0^2-2g\left(z-z_0\right)$
Frictionless rolling (rigid body, $\beta =I/m{R}^2$ )	$\frac{-g\text{d}z/\text{d}s}{1+\beta }$	$V_0^2-\frac{2g\left(z-z_0\right)}{1+\beta }$
Coulomb friction ($\mu$ = friction coeff.)	$-g\frac{\text{d}z}{\text{d}s}-\mu \mathrm{sgn}\left(V\right)\mathrm{sgn}\left(N\right)\times \left(g\frac{\text{d}x}{\text{d}s}+\kappa V^2\right)$	${\text{e}}^{-2\mu \Theta \left(s\right)}\left[V_0^2+2{\displaystyle {\int }_{s_0}^s{\text{e}}^{2\mu \Theta \left(s^{\prime }\right)}{A}_0\left(s^{\prime }\right)\text{d}{s}^{\prime }}\right]$
Viscous (quadratic) drag ($b/m$ = drag/mass ratio)	$-g\frac{\text{d}z}{\text{d}s}-\frac{b}{m}{V}^2$	First-order linear ODE in $V^2$ ; integrating factor ${\text{e}}^{2\left(b/m\right)\left(s-s_0\right)}$ :${\text{e}}^{-2\frac{b}{m}\left(s-s_0\right)}\left[V_0^2+2{\displaystyle {\int }_{s_0}^s{\text{e}}^{2\frac{b}{m}\left(s^{\prime }-s_0\right)}{\mathcal{K}}_g\left(s^{\prime }\right)\text{d}{s}^{\prime }}\right]$
Droplet (surface tension) ($\gamma$ = surface tension, ${\ell }_c$ = contact length)	$-g \frac{\text{d}z}{\text{d}s}-\frac{\gamma {\ell }_c\kappa }{m}\text{sgn}\left(V\right)$	Numerical; closed form only forconstant-curvature curves
Rotating frame (Ω = angular velocity)	$-g\frac{\text{d}z}{\text{d}s}+{\Omega }^{2}x\frac{\text{d}x}{\text{d}s}$	$V_0^2-2g\left(z-z_0\right)+{\Omega }^2\left(x^2-x_0^2\right)$

(A4) Existence of turning points.

Proposition 2 locates a turning point as a root of $z\left(s_{\star }\right)=z_0+V_0^2/\left(2g\right)$ . Such a root exists within the curve’s domain if and only if $z\left(\lambda \right)$ attains the value $z_0+V_0^2/\left(2g\right)$ for some $\lambda$ in the admissible range; this holds, for example, whenever $z\left(\lambda \right)$ is continuous and the curve is bounded above by at least this height (true for all eight families of Table 3 on a sufficiently large domain). If no such root exists—e.g. an unbounded monotonic ramp with $V_0^2/\left(2g\right)$ exceeding the total rise available—the particle never turns back, $V\left(s\right)>0$ throughout, and Pass 2’s integration domain $\left[s_0,\infty \right)$ must be treated as open-ended rather than bounded by $s_{\star }$ ; Pass 3 tiling is then inapplicable since there is no periodic motion to tile. For the dissipative kernels of Section 3, the analogous turning point must instead be obtained as the first zero of $V^2\left(s\right)$ from the kernel ODE (e.g. (23)) rather than from (14), as already remarked after Proposition 2.

Table 3. Dimensional reduction: arc-length formulation specialised to eight curve families. ${\phi }^{\prime }$ denotes $\text{d}\phi /\text{d}x$ ; $\dot{x}=\text{d}x/\text{d}\lambda$; $r^{\prime }=\text{d}r/\text{d}\theta$ ; $G=F_{xx}{F}_z^2-2F_{xz}{F}_x{F}_z+F_{zz}{F}_x^2$ for implicit curves.

Geometry	$\lambda$	$\text{d}s/\text{d}\lambda$	$z\left(\lambda \right)$	$\kappa \left(\lambda \right)$
Cartesian graph $z=\phi \left(x\right)$	$x$	$\sqrt{1+{{\phi }^{\prime }}^2}$	$\phi \left(x\right)$	${\phi }^{″}/{\left(1+{{\phi }^{\prime }}^2\right)}^{3/2}$
Circle/Pendulum $r=L$ ($\theta$ from downward vertical; $z=0$ at pivot, min at $\theta =0$ )	$\theta$	$L$	$-L\cos \theta$	$1/L$
Ellipse $\frac{x^2}{a^2}+\frac{z^2}{b^2}=1$	$\theta$ (eccentric)	$\sqrt{a^2{\sin }^2\theta +b^2{\cos }^2\theta }$	$b\sin \theta$	$ab/{\left(a^2{\sin }^2\theta +b^2{\cos }^2\theta \right)}^{3/2}$
Parabola $z=a{x}^2$	$x$	$\sqrt{1+4a^2{x}^2}$	$a{x}^2$	$2a/{\left(1+4a^2{x}^2\right)}^{3/2}$
Log. Spiral $r=r_0{\text{e}}^{b\theta }$	$\theta$	$r\sqrt{1+b^2}$	$r\sin \theta$	$1/\left(r\sqrt{1+b^2}\right)$
Catenary $z=c\cosh \left(x/c\right)$	$x$	$\text{cosh}\left(x/c\right)$	$c\cosh \left(x/c\right)$	$1/\left(c{\cosh }^2\left(x/c\right)\right)$
Cycloid $x=R\left(\theta -\sin \theta \right)$ , $z=R\left(1-\cos \theta \right)$	$\theta$	$2R\sin \left(\theta /2\right)$	$R\left(1-\cos \theta \right)$	$-1/\left(4R\sin \left(\theta /2\right)\right)$
Implicit $F\left(x,z\right)=0$	$s$	1	$z\left(s\right)$	$G/{\left|\nabla F\right|}^3$

2.4. Force Decomposition

The gravitational force per unit mass is $g=\left(0,-g\right)$ . Its projections onto the Frenet-Serret frame are:

$g_{\parallel }\left(s\right)=g\cdot \widehat{t}=-g\frac{\text{d}z}{\text{d}s},$(9)

$g_{\perp }\left(s\right)=g\cdot \widehat{n}=-g\frac{\text{d}x}{\text{d}s}.$(10)

Equation (9) embodies the key simplification: the tangential driving term is simply $-g\text{d}z/\text{d}s$ , the rate of height change with arc length—two symbols replacing the Cartesian $-g{\phi }^{\prime }/\left(1+{{\phi }^{\prime }}^2\right)$ .

2.5. Equation of Motion and the Spatial Speed Field

Newton’s second law projected onto $\widehat{t}$ (frictionless case) gives

$m\frac{{\text{d}}^{2}s}{\text{d}{t}^2}=-mg\frac{\text{d}z}{\text{d}s}.$ (11)

Let $V\left(s\right)=\text{d}s/\text{d}t$ . Then ${\text{d}}^{2}s/\text{d}{t}^2=V\text{d}V/\text{d}s$ , so (11) becomes

$V\frac{\text{d}V}{\text{d}s}=-g\frac{\text{d}z}{\text{d}s}.$ (12)

Integrating from $s_0$ to $s$ :

Theorem 1 (Spatial Speed Field). For frictionless, holonomic constrained motion on a smooth planar curve under gravity, the arc-length speed satisfies

$\boxed{V^2\left(s\right)=V_0^2-2g\left[z\left(s\right)-z_0\right],}$ (13)

where $z_0=z\left(s_0\right)$ is the initial height.

Proof. Direct integration of (12) from $s_0$ to $s$ gives $\frac{1}{2}{V}^2\left(s\right)-\frac{1}{2}{V}_0^2=-g\left[z\left(s\right)-z_0\right]$ .

Remark. Equation (13) requires only the height function $z\left(s\right)$ and involves no arc-length quadrature. For a curve parametrized by a generic parameter $\lambda$ , one simply substitutes $z\left(\lambda \right)$ directly; the arc-length element $\text{d}s/\text{d}\lambda$ enters only in Pass 2.

Proposition 2 (Turning Points). A turning point $s_{\star }$ satisfies

$z\left(s_{\star }\right)=z_0+\frac{V_0^2}{2g},$ (14)

independent of parametrization and dependent only on height.

Remark. Equation (14) is a pure Pass-1 result: it requires only the height function $z\left(\lambda \right)$ and the initial conditions, involving no quadrature. In practice it serves two purposes. First, it provides a geometry-based upper bound on the reachable arc, allowing Pass-2 to be restricted to a finite integration domain $\left[s_0,s_{\star }\right]$ without searching for the zero of $V^2$ during time integration. Second, for the frictionless case it furnishes an independent check on the Pass-2 turning-point location computed by root-finding on the speed field. Note that (14) applies only when gravity is the sole force; with dissipative forces such as Coulomb friction the turning point must be determined from the appropriate kernel ODE (Section 3), as in Example 2.

2.6. Normal Contact Force

Newton’s second law projected onto $\widehat{n}$ gives

Theorem 3 (Normal Contact Force).

$\boxed{N\left(s\right)=m\left[g\frac{\text{d}x}{\text{d}s}+\kappa \left(s\right)V^2\left(s\right)\right].}$ (15)

Proof. Newton’s second law projected onto $\widehat{n}$ gives the centripetal balance $m\kappa \left(s\right) V^2\left(s\right)=N+m g\cdot \widehat{n}$ . Substituting $g\cdot \widehat{n}=-g\text{d}x/\text{d}s$ from (10) and rearranging yields (15).

The first term $mg\text{d}x/\text{d}s$ is the static weight component perpendicular to the slope; the second term $m\kappa V^2$ is the dynamic centripetal contribution. Both are computed from Pass-1 quantities alone. For a bilateral constraint, $N$ may be positive (inward) or negative (outward); the particle remains on the track in both cases.

Corollary 1 (Force Sign and Bilateral Constraint). For a unilateral constraint, contact is lost when $N\left(s\right)=0$ :

$g\frac{\text{d}x}{\text{d}s}+\kappa \left(s\right)V^2\left(s\right)=0.$ (16)

Substituting (13) transforms (16) into an algebraic equation in $s$ solvable from Pass-1 data. For bilateral constraints, (16) locates the transition between inward and outward rail loading.

2.7. The Three-Pass Solution Pipeline

2.7.1. Pass 1—Spatial Quantities (Exact)

Given the curve $r\left(\lambda \right)$ and initial conditions $\left(s_0,V_0\right)$ :

1) Compute $\kappa \left(\lambda \right)$ , $\text{d}x/\text{d}s$ , $\text{d}z/\text{d}s$ from (8) and (6).

2) Evaluate $V^2\left(\lambda \right)$ from (13)—exact, no quadrature.

3) Evaluate $N\left(\lambda \right)$ from (15)—exact.

4) Solve (14) for turning point(s)—algebraic.

All Pass-1 results are closed-form elementary functions whenever $z\left(\lambda \right)$ is.

2.7.2. Pass 2—Temporal Quadrature

$t\left(s\right)={\displaystyle {\int }_{s_0}^s\frac{\text{d}{s}^{\prime }}{V\left(s^{\prime }\right)}}={\displaystyle {\int }_{{\lambda }_0}^{\lambda }\frac{\left|\text{d}r/\text{d}{\lambda }^{\prime }\right|\text{d}{\lambda }^{\prime }}{\sqrt{V_0^2-2g\left[z\left({\lambda }^{\prime }\right)-z_0\right]}}}.$ (17)

The integral (17) is non-elementary in general—for a circular arc it reduces to the elliptic integral $K\left(k\right)$ , which is irreducibly non-elementary (Abel’s result for elliptic integrals; the general algebraic case is covered by Liouville’s theorem [13]). It is evaluated by adaptive Gauss-Kronrod quadrature [14] to machine precision. This choice is well-suited to the integrand structure: the integrand $\left|\text{d}r/\text{d}{\lambda }^{\prime }\right|/\sqrt{V_0^2-2g\left[z\left({\lambda }^{\prime }\right)-z_0\right]}$ is smooth on the interior of each half-period arc but develops an integrable algebraic singularity of the form ${\left(s_{\star }-s\right)}^{-1/2}$ near turning points; Gauss–Kronrod rules with adaptive subdivision concentrate nodes near such singularities automatically, achieving machine precision with a modest number of function evaluations without requiring a change of variable or special endpoint treatment. The half-period $\tau =t\left(s_{\star }\right)$ is computed once and stored; this formula gives the true half-period only when the initial point $s_0$ is itself a turning point ($V_0=0$ ). For $V_0>0$ , the half-period must instead be computed as the time-of-flight between the two bounding turning points of the orbit.

2.7.3. Pass 3—Inversion and Tiling

Since $t\left(s\right)$ is strictly monotone on each arc between turning points, the inverse $s\left(t\right)$ is computed by interpolation on the precomputed Pass-2 table. For conservative (periodic) systems, global trajectories for large $t$ are recovered by exact tiling:

$s\left(t\right)=s\left(t\mathrm{mod}\tau \right)$ , with arc direction reversed after each half-period.(18)

The tiling phase error in (18) is bounded by $N{\epsilon }_{\tau }$ , where ${\epsilon }_{\tau }$ is the quadrature precision of $\tau$ ; tiling introduces no additional quadrature at large time (see Section 6.6 for the full error bound). For dissipative systems the motion is non-periodic, tiling is inapplicable, and $s\left(t\right)$ is obtained by direct interpolation on the Pass-2 table alone. This contrasts with direct ODE integration, where phase error accumulates as $O\left(h^p T_{\text{final}}\right)$ for a $p$ -th order method with step size $h$ .

3. Speed-Field Kernels

The tangential equation of motion (11) can be written in the unified form

$V \frac{\text{d}V}{\text{d}s}=\mathcal{A}\left(s;V\right),$ (19)

where the speed-field kernel $\mathcal{A}$ encodes the physical force model. Structurally, $\mathcal{A}$ acts as an influence kernel for the speed-squared field: it specifies the contribution of each arc-length element $\text{d}s$ to the total speed change, in exact analogy with the discrete convolution kernel of (1). When $\mathcal{A}$ is independent of $V$ , (19) integrates to a closed-form $V^2\left(s\right)$ . When $\mathcal{A}$ depends linearly on $V^2$ , it becomes a first-order linear ODE in $V^2$ with an integrating-factor solution.

The kernel structure is the key to extensibility: any force law that can be expressed as a function of arc-length position (and possibly $V^2$ ) yields a modified kernel $\mathcal{A}$ , while the downstream Pass-2 and Pass-3 steps remain unchanged. Table 2 summarises the kernels and their solutions for six models.

Note on Coulomb row: $A_0\left(s\right)=-g\text{d}z/\text{d}s-\mu \text{sign}\left(V\right)\text{sign}\left(N\right)g\text{d}x/\text{d}s$ is the part of the kernel that does not multiply $V^2$ , and $\text{Θ}\left(s\right)={\displaystyle {\int }_{s_0}^s\kappa \text{d}{s}^{\prime }}$ is the accumulated turning angle from the initial position. The closed-form solution in the $V^2$ column is piecewise: it holds on each sub-arc $\left[s_{*},s_{**}\right]$ over which $\text{sign}\left(N\right)$ does not change and the direction of motion $\text{sign}\left(V\right)$ is constant (the case $\text{sign}\left(V\right)=\text{sign}\left(N\right)=+1$ underlying (24)); if $\text{sign}\left(V\right)\text{sign}\left(N\right)=-1$ on a sub-arc, the sign in the exponent of the integrating factor flips accordingly. Here $V_{*}^2$ is the speed-squared at the left endpoint $s_{*}$ of that sub-arc. At each crossing where $g\text{d}x/\text{d}s+\kappa V^2=0$ (i.e. $N=0$ ), the sign of $N$ must be re-evaluated from (23) and the integration restarted with the current state as the new initial condition. The general multi-sign-change algorithm is described in Section 3.2 (Equations (22)-(23)).

The rotating-frame kernel in the last row illustrates how non-gravitational potential fields are accommodated: the term ${\Omega }^{2}x\text{d}x/\text{d}s$ is simply the centrifugal potential gradient projected onto the tangent, and it contributes additively to the kernel $\mathcal{A}$ alongside gravity. Any conservative force whose potential $U\left(x,z\right)$ is known can be incorporated by replacing the gravitational term $g\text{d}z/\text{d}s$ with $-\nabla U\cdot \widehat{t}=-\left(U_x\text{d}x/\text{d}s+U_z\text{d}z/\text{d}s\right)$, yielding $V^2\left(s\right)=V_0^2+2\left[U\left(x_0,z_0\right)-U\left(x\left(s\right),z\left(s\right)\right)\right]$ as the generalized Pass-1 result. The quadrature pipeline of Pass 2 and tiling of Pass 3 then proceed unchanged.

3.1. Rolling with Inertia

For a rigid body rolling without slipping, the moment of inertia $I$ modifies the effective mass through the Lagrange-d’Alembert constraint:

$\left(m+I/R^2\right)V\frac{\text{d}V}{\text{d}s}=-mg\frac{\text{d}z}{\text{d}s},$ (20)

so that $V^2\left(s\right)=V_0^2-2g\left(z-z_0\right)/\left(1+\beta \right)$ with $\beta =I/m{R}^2$ , where $R$ is the radius of the rolling body itself, not the radius of curvature of the track. For a solid sphere $\beta =2/5$ ; for a solid cylinder $\beta =1/2$ ; for a thin ring $\beta =1$ .

3.2. Coulomb Friction

With kinetic friction force $f=\mu \left|N\right|$ opposing motion, (11) becomes

$V\frac{\text{d}V}{\text{d}s}=-g\frac{\text{d}z}{\text{d}s}-\mu \text{sign}\left(V\right)\frac{\left|N\left(s\right)\right|}{m}.$ (21)

From (15), $N\left(s\right)=m\left[g\text{d}x/\text{d}s+\kappa \left(s\right)V^2\left(s\right)\right]$ , which may be positive or negative for a bilateral constraint. Since $\left|N\right|/m=\left|g\text{d}x/\text{d}s+\kappa \left(s\right)V^2\right|$ , the friction term is most transparently written without introducing $\text{sign}\left(N\right)$ :

$\frac{\text{d}\left(V^2\right)}{\text{d}s}=-2g\frac{\text{d}z}{\text{d}s}-2\mu \text{sign}\left(V\right)\left|g\frac{\text{d}x}{\text{d}s}+\kappa \left(s\right)V^2\right|.$ (22)

Equation (22) is valid for all signs of $N$ and $V$ without restriction. To convert it to a linear ODE in $V^2$ one writes $\left|g\text{d}x/\text{d}s+\kappa V^2\right|=\text{sign}\left(N\right)\left(g\text{d}x/\text{d}s+\kappa V^2\right)$ , which requires $\text{sign}\left(N\right)\left(g\text{d}x/\text{d}s+\kappa V^2\right)\ge 0$ . This identity holds whenever the contact force computed from Pass 1 correctly predicts the sign of $N$ , i.e. whenever $g\text{d}x/\text{d}s+\kappa V^2$ and $N$ share the same sign. If the particle approaches liftoff ($N\to 0$ ) or if $\kappa V^2$ becomes large enough to change the bracket’s sign, $\text{sign}\left(N\right)$ must be re-evaluated from the current state before advancing. Subject to this condition, substituting $|\cdot |=\text{sign}\left(N\right)(\cdot )$ and multiplying through by $2V$ gives

$\begin{array}{l}\frac{\text{d}\left(V^2\right)}{\text{d}s}=-2g\frac{\text{d}z}{\text{d}s}-2\mu \text{sign}\left(V\right)\text{sign}\left(N\right)\left[g\frac{\text{d}x}{\text{d}s}+\kappa \left(s\right)V^2\right],\\ \text{sign}\left(N\right)\left[g\frac{\text{d}x}{\text{d}s}+\kappa V^2\right]\ge 0.\end{array}$ (23)

This is a first-order linear ODE in $V^2\left(s\right)$ . In the practically common case where contact is inward ($N>0$ , so $\text{sign}\left(N\right)=+1$ ) and motion is in the positive-$s$ direction ($\text{sign}\left(V\right)=+1$ ), (23) reduces to

$\frac{\text{d}\left(V^2\right)}{\text{d}s}=-2g\left(\frac{\text{d}z}{\text{d}s}+\mu \frac{\text{d}x}{\text{d}s}\right)-2\mu \kappa \left(s\right)V^2,$ (24)

which has the integrating-factor solution shown in Table 2, with integrating factor $\exp \left(2\mu \Theta \left(s\right)\right)$ where $\Theta \left(s\right)={\displaystyle {\int }_{s_0}^s\kappa \text{d}{s}^{\prime }}$ is the total turning angle (equal to 2π for a closed convex curve, by the Gauss–Bonnet theorem). In Example 2, $\text{sign}\left(N\right)$ changes at $x\approx -0.1052$ m; the pipeline monitors $g\text{d}x/\text{d}s+\kappa V^2$ at each quadrature node and updates $\text{sign}\left(N\right)$ accordingly, so (23) remains valid throughout.

Practical Procedure for Repeated Sign Changes and Liftoff

Equation (23) is piecewise-valid: it holds on each maximal sub-arc over which $\text{sign}\left(N\right)$ is constant, and must be restarted whenever that sign changes. In a track with multiple inflections or oscillatory normal loading, $N\left(s\right)$ may cross zero more than once, so the following stepwise rule is applied uniformly (this is the procedure formalized in the Pass-2 sign-monitoring/restart steps of Algorithm 1):

Algorithm 1. The GQM three-pass pipeline.

1) Initialize. At the current sub-arc’s left endpoint $s_{*}$ , evaluate $\sigma \equiv \text{sign}\left[g\text{d}x/\text{d}s+\kappa \left(s_{*}\right)V^2\left(s_{*}\right)\right]$ and use this as $\text{sign}\left(N\right)$ in (23) for the remainder of the sub-arc.

2) Monitor. While integrating (23) forward (by quadrature or by the equivalent linear-ODE closed form), evaluate the bracketed quantity $g\text{d}x/\text{d}s+\kappa \left(s\right)V^2\left(s\right)$ at each quadrature node (or at each adaptive step).

3) Detect a sign change. If the bracketed quantity’s sign disagrees with $\sigma$ at some node, a crossing $N\left(s_c\right)=0$ has occurred between the previous and current node. Bracket $s_c$ between these two nodes and locate it by root-finding (bisection or Brent’s method) on $g\text{d}x/\text{d}s+\kappa \left(s\right)V^2\left(s\right)=0$ , using the already-computed $V^2\left(s\right)$ on that sub-arc.

4) Restart. Terminate the current sub-arc integration at $s_c$ . Set the new initial condition $\left(s_0,V_0^2\right)\leftarrow \left(s_c,V^2\left(s_c\right)\right)$ , recompute $\sigma$ from the (now opposite) sign of the bracketed quantity just past $s_c$ , and resume integration of (23) with the updated $\sigma$ as the new $\text{sign}\left(N\right)$ . Repeat from Step (2).

5) Liftoff (unilateral constraints only). If the constraint is unilateral and the bracketed quantity approaches zero from the positive (inward-contact) side without the centripetal term $\kappa V^2$ being able to restore $N>0$ immediately afterward—i.e. the particle would require $N<0$ to remain on the track—then contact is lost at $s_c$ rather than merely changing sign. In this case the friction force vanishes identically beyond $s_c$ (there is no normal load to generate it), the particle leaves the curve, and its subsequent motion is governed by free-flight (projectile) kinematics under gravity alone, decoupled from the curve’s geometry. This regime is outside the scope of the constrained-motion pipeline of Section 2.7 and is not pursued further here (see Section 6.7). For a bilateral constraint, by contrast, no such hand-off is needed: $N$ is simply allowed to become negative (outward loading), and Step (4) applies as stated.

For practical implementation, Step (3)’s bracket search is most robust when performed on the dense/continuous output of the Pass-2 quadrature (or an auxiliary fine grid) rather than relying solely on the adaptive quadrature’s own node spacing, since Gauss-Kronrod nodes are optimized for integral accuracy rather than for resolving zero crossings of $N\left(s\right)$ to high precision.

4. Dimensional Reduction Table

The arc-length formulation specialises to each coordinate system by substituting the appropriate expressions for $z\left(\lambda \right)$ , $\text{d}s/\text{d}\lambda$ , and $\kappa \left(\lambda \right)$ . Table 3 collects these for eight standard curve families. In every case Pass 1 is exact: $V^2\left(\lambda \right)=V_0^2-2g\left(z\left(\lambda \right)-z_0\right)$ requires only $z\left(\lambda \right)$ . The arc-length element $\text{d}s/\text{d}\lambda$ appears only in the Pass-2 time integral (17).

5. Numerical Examples

The GQM pipeline is evaluated on two representative test cases: a simple pendulum loop and a cubic curve with an inflection point and Coulomb friction. For each case, GQM results are compared against conventional time-stepping schemes. Table 4 summarizes key accuracy metrics for the long-duration pendulum simulation.

5.1. Example 1: Simple Pendulum Loop

We first benchmark GQM against a fixed-step fourth-order Runge-Kutta scheme (RK4, $h=0.01$ s) for a simple pendulum over a long simulation window $T_{\text{max}}=1000\text{s}$ . The pendulum has length $L=1\text{m}$ and is released from rest ($V_0=0$ ) at an initial angle ${\theta }_0=-\pi /2+1.2\approx -0.371$ rad, which is 1.2 rad above the equilibrium position ${\theta }_{\text{eq}}=-\pi /2$ . The pendulum angle $\theta$ is measured from the horizontal, so the pivot is at the origin, the lowest position corresponds to ${\theta }_{\text{eq}}=-\pi /2$ , and the height function is $z\left(\theta \right)=L\sin \theta$ (not $-L\cos \theta$ as in the

Table 4. Numerical precision metric comparing GQM against explicit RK4 (simple pendulum, $L=1\text{m}$ , $V_0=0$ , $T_{\text{max}}=1000\text{s}$ , $h=0.01\text{s}$ ). This is the single metric directly computed and reported by the supplementary script GQM_Plots.py; see the discussion following the table for the analytical bound governing GQM’s expected drift on metrics not computed here.

Performance metric	Explicit RK4 ($h=0.01\text{s}$ ), measured vs. GQM reference	Proposed GQM(self-consistency)
Angular phase drift $\text{Δ}\theta$ (at $t={10}^3\text{s}$ )	1.0 × 10−4 rad	not an independent angular measurement*

*GQM’s angular drift is not measured against an independent reference in this experiment (GQM is itself the reference against which RK4 is compared); rather, GQM reconstructs $\theta \left(t\right)$ analytically from the closed-form speed field (25) via exact half-period tiling (18), so its self-consistency error is bounded by the time-domain quantity $N{\epsilon }_{\tau }+{\epsilon }_{\text{interp}}\left(\delta t\right)\approx 5\times {10}^{-13}\text{s}$ (a time, not an angle; see (32) and Section 6.6 for the derivation and its units), evaluated at $N\approx 454$ periods with ${\epsilon }_{\tau }\approx {10}^{-14}\text{-}{10}^{-15}\text{s}$ and ${\epsilon }_{\text{interp}}\approx {10}^{-14}\text{s}$ . This is not an independently measured angular residual.

angle-from-vertical convention of Table 3). The two conventions describe the same physical circle; here $z\left(\theta \right)=L\sin \theta$ is used throughout so that equilibrium at $\theta =-\pi /2$ gives $z_{\text{eq}}=-L$ , consistent with the figure. This value is used consistently in the text, figure, and all numerical results; the full simulation parameters are also encoded in the supplementary script GQM_Plots.py. For reproducibility, the Pass-1 speed field and Pass-2 integrand are stated explicitly:

$V^2\left(\theta \right)=V_0^2-2g\left[L\sin \theta -L\sin {\theta }_0\right]=-2gL\left(\sin \theta -\sin {\theta }_0\right),$ (25)

and the Pass-2 time-of-flight integral is

$t\left(\theta \right)={\displaystyle {\int }_{{\theta }_0}^{\theta }\frac{L\text{d}{\theta }^{\prime }}{\sqrt{-2gL\left(\sin {\theta }^{\prime }-\sin {\theta }_0\right)}}},$ (26)

with the half-period $\tau$ computed as $t\left({\theta }_{\star }\right)$ where ${\theta }_{\star }$ satisfies $V^2\left({\theta }_{\star }\right)=0$ , i.e. $\sin {\theta }_{\star }=\sin {\theta }_0+V_0^2/\left(2gL\right)=\sin {\theta }_0$ . The equation $\sin \theta =\sin {\theta }_0$ admits two solutions: the trivial root $\theta ={\theta }_0$ and the physical (reachable) turning point ${\theta }_{\star }=-\pi -{\theta }_0\approx -2.771\text{rad}$ , which is the point symmetric to ${\theta }_0$ about the bottom ${\theta }_{\text{eq}}=-\pi /2$ . The arc-length element for the circle of radius $L$ is $\text{d}s=L\text{d}\theta$ , which has been substituted directly into (26). All numerical results in Table 4 and Figure 1 were obtained using equations (25)-(26).

As shown in Figure 1 and Table 4, conventional RK4 integration exhibits gradual phase drift and amplitude decay due to truncation error. At $t=1000\text{s}$ , the RK4 solution differs from the GQM reference by $\text{Δ}\theta =0.0001\text{rad}$ . GQM’s own self-consistency error is bounded by the time-domain quantity $N{\epsilon }_{\tau }+{\epsilon }_{\text{interp}}\left(\delta t\right)\approx 5\times {10}^{-13}\text{s}$ (see (32)); converting this to an angular-equivalent bound via $\text{Δ}\theta ~\dot{\theta }\left(N{\epsilon }_{\tau }+{\epsilon }_{\text{interp}}\right)$, with $\dot{\theta }=V/L$ of order 1 rad/s near the bottom of the swing, gives an angular drift on the order of 10⁻¹² - 10⁻¹³ rad—many orders of magnitude below the 10⁻⁴ rad RK4 drift. Because GQM reconstructs the motion

Figure 1. Benchmark verification for a simple pendulum loop ($T_{\text{max}}=1000\text{s}$ ; release angle ${\theta }_0=-\pi /2+1.2\approx -0.371\text{rad}$ , corresponding to a displacement of 1.2 rad above the equilibrium position ${\theta }_{\text{eq}}=-\pi /2$ ). The top-left panel shows phase alignment over the initial swings. The top-right panel compares the GQM and RK4 solutions at $t=1000\text{s}$ and highlights the accumulated phase difference $\text{Δ}\theta =0.0001\text{rad}$ . The bottom panels show the instantaneous relative energy variation and the cumulative energy error, demonstrating that GQM maintains energy to machine precision over the full simulation.

using exact half-period tiling based on high-accuracy spatial quadrature, it reduces phase and energy drift to this level, maintaining near-machine-precision agreement over the full time horizon.

5.2. Example 2: Cubic Curve with Inflection Point and Coulomb Friction

5.2.1. Physical Setup

A unit-mass particle slides along the cubic track $z=x^3$ from an initial point $M_0=\left(-1,-1\right)$ with starting speed $V_0=6.2088\text{m}\cdot {\text{s}}^{-1}$ and Coulomb friction coefficient $\mu =0.1$ . The track passes through the origin, where the signed curvature $\kappa \left(s\right)$ changes sign: the segment in the third quadrant ($x<0$ ) is convex, while the segment in the first quadrant ($x>0$ ) is concave. This curvature sign reversal constitutes an inflection point and represents a scenario where conventional Frenet-frame implementations produce a 180˚ flip in the normal direction, introducing artificial discontinuities in the computed normal reaction force.

In GQM the normal $\widehat{n}\left(s\right)$ is defined by a fixed global rotation of the tangent (cf. (7)), so it remains continuous and well-defined through the inflection point. The signed curvature $\kappa \left(s\right)$ transitions smoothly through zero at the origin without any discontinuity in the frame or loss of numerical stability.

5.2.2. Pass-1 Spatial Quantities (Exact)

Using the Cartesian-graph row of Table 3 with $\phi \left(x\right)=x^3$ :

$\frac{\text{d}s}{\text{d}x}=\sqrt{1+9x^4},$ (27)

$\kappa \left(x\right)=\frac{6x}{{\left(1+9x^4\right)}^{3/2}}.$ (28)

The speed field under Coulomb friction is governed by the linear ODE (23) with $\text{d}z/\text{d}s=3x^2/{\left(1+9x^4\right)}^{1/2}$ and $\text{d}x/\text{d}s=1/{\left(1+9x^4\right)}^{1/2}$ . At the inflection point $x=0$ one has $\kappa =0$ , so the normal force reduces to $N/m=g{\text{d}x/\text{d}s|}_{x=0}=g$ (since ${\text{d}x/\text{d}s|}_{x=0}=1$ ). The kernel (23) then gives

$\begin{array}{c}{\frac{\text{d}\left(V^2\right)}{\text{d}s}|}_{x=0}=-2g{\frac{\text{d}z}{\text{d}s}|}_{x=0}-2\mu \text{sign}\left(V\right)\text{sign}\left(N\right)g{\frac{\text{d}x}{\text{d}s}|}_{x=0}\\ =-2\mu \text{sign}\left(V\right)\text{sign}\left(N\right)g,\end{array}$ (29)

using ${\text{d}z/\text{d}s|}_{x=0}=0$ . This is finite and continuous—confirming that no special treatment is needed at the inflection point.

5.2.3. Turning Point and Normal-Force Sign Change

The particle decelerates due to both the rising elevation and friction, coming to rest at the dissipative turning point where $V\left(s_{\star }\right)=0$ . This turning point is located at

$\begin{array}{l}{M}_{\star }=\left(x_{\star },z_{\star }\right)\approx \left(0.8367,0.5857\right)\text{m},\\ s_{\star }-s_0={\displaystyle {\int }_{-1}^{x_{\star }}\sqrt{1+9x^4}\text{d}x}\approx 2.6502\text{m},\end{array}$ (30)

computed to machine precision using a single adaptive Gauss–Kronrod evaluation, where $s_0$ denotes the arc-length coordinate at $M_0=\left(-1,-1\right)$ . The normal reaction $N\left(s\right)$ changes sign at $x\approx -0.1052\text{m}$ , $z\approx -0.0012\text{m}$ , just before the inflection point; for a unilateral constraint this is the location at which contact would be lost.

5.2.4. Results and Figure

Figure 2 shows the geometric path and the continuous signed curvature profile $\kappa \left(s\right)$ along the cubic track. The left panel displays the particle trajectory from $M_0\left(-1,-1\right)$ through the inflection at the origin to the turning point $M_{\star }$ , colour-coded by the sign of the normal reaction $N\left(s\right)$ . The right panel shows $\kappa \left(s\right)$ passing smoothly through zero at $x=0$ with no discontinuity, providing a stable and continuous input to the dissipative kernel throughout the motion. Table 5 summarises the corresponding Pass-1 numerical results.

Figure 2. Example 2: Cubic track $z=x^3$ with inflection point and Coulomb friction ($\mu =0.1$ , $V_0=6.2088\text{m}\cdot {\text{s}}^{-1}$ , $M_0=\left(-1,-1\right)$ , $m=1\text{kg}$ , $g=9.81\text{m}\cdot {\text{s}}^{-2}$ ). Left panel—geometric path from $M_0\left(-1,-1\right)$ through the origin to the dissipative turning point $M_{\star }\approx \left(0.837,0.586\right)$ ; colour coding indicates the sign of the normal reaction $N\left(s\right)$ . Right panel—continuous signed curvature $\kappa \left(s\right)$ along the track, passing smoothly through zero at $x=0$ (the inflection point) and providing a stable input to the dissipative kernel with no frame-flip artefacts.

Table 5. Example 2 (Cubic curve with inflection)—numerical results. All Pass-1 quantities are exact closed-form results from (13) and (15) with $\mu =0.1$ , $V_0=6.2088\text{m}\cdot {\text{s}}^{-1}$ , $m=1\text{kg}$ , $g=9.81\text{m}\cdot {\text{s}}^{-2}$ .

Metric	Value	Unit
Turning point position $M_{\star }=\left(x_{\star },z_{\star }\right)$	(0.8367, 0.5857)	m
Arc length $M_0\to M_{\star }$ : ${\displaystyle {\int }_{x_0}^{x_{\star }}\sqrt{1+9x^4}\text{d}x}$	2.6502	m
$N=0$ crossing position $x_{N=0}$	−0.1052	m
$\kappa$ at inflection ($x=0$ )	0 (exact)	m−1
Frame continuity through inflection	Continuous (no flip)	-

5.2.5. Independent Verification by Direct Time-Domain Integration

The GQM results above are obtained by solving the spatial kernel ODE (23) in the arc-length (equivalently, $x$ ) domain, with the sign-monitoring and restart procedure of Section 3.2. As an independent check, we additionally integrate the same physical model—Newton’s second law (21) with the friction law $f=\mu \left|N\right|$ —directly in the time domain, using the state vector $\left(x\left(t\right),V\left(t\right)\right)$ and a standard adaptive Runge-Kutta scheme (Dormand-Prince RK45), entirely independent of the spatial quadrature pipeline. The two formulations share only the underlying physics, not the numerical method: the spatial integration of Section 3.2 advances $V^2$ as a function of $x$ via quadrature with explicit sign-of-$N$ bookkeeping, whereas the time-domain check advances $\left(x,V\right)$ as a function of $t$ via a general-purpose ODE solver with event detection for $V=0$ and a posteriori root-finding for $N=0$ .

Table 6 compares the two computations. The turning-point location and the $N=0$ crossing agree to within 10⁻⁶ m (limited by the time-domain integrator’s step-size tolerance, $\text{Δ}{t}_{\text{max}}={10}^{-4}\text{s}$ ); this agreement confirms that the spatial pipeline’s sign-monitoring and restart logic (Section 3.2, Algorithm 1) correctly reproduces the physical model rather than an artefact of the spatial-quadrature implementation.

Table 6. Example 2 cross-check: GQM spatial pipeline vs. independent direct time-domain integration (Dormand-Prince RK45, adaptive step, event-based turning-point and $N=0$ detection) of the same Coulomb-friction model (21).

Quantity	GQM (spatial pipeline)	Time-domain RK45	Difference
Turning point $x_{\star }$ (m)	0.83669	0.83669	1.2 × 10−6
Turning point $z_{\star }$ (m)	0.58573	0.58573	2.5 × 10−6
Arc length to $M_{\star }$ (m)	2.65020	2.65020*	2.7 × 10−6
$N=0$ crossing $x_{N=0}$ (m)	−0.10523	−0.10523	3.1 × 10−10
Time of flight to $M_{\star }$ (s)	-†	0.8203	-

*Arc length along the time-domain trajectory is recovered by evaluating (17)’s arc-length integral ${\displaystyle \int \sqrt{1+9x^4}\text{d}x}$ on the RK45-computed $x\left(t\right)$ path. ^†The spatial pipeline does not require the time-of-flight value to obtain $M_{\star }$ or the $N=0$ crossing (these are recovered purely from Pass-1 and the quadrature in $x$ ); the time-domain figure is reported here only as an additional consistency datum, since it falls out naturally from the time-domain solve.

6. Discussion

6.1. Advantages of the Arc-Length Formulation

The primary advantage is the elimination of the $\sqrt{1+{{\phi }^{\prime }}^2}$ factor, which diverges at vertical tangents and prevents the Cartesian formulation from treating the upper arc of a circle, ellipse sides, or any polar curve where $\left|{\phi }^{\prime }\right|\to \infty$ . The arc-length formulation is coordinate-free and applies to any rectifiable plane curve.

A secondary advantage is algebraic transparency. The speed field (13) states in three symbols that squared speed equals its initial value minus twice the height gain. The normal force (15) is similarly compact: $g\text{d}x/\text{d}s$ is the fraction of gravitational weight supported by the surface, and $\kappa V^2$ is the centripetal demand. Both derive from Newton’s second law—energy conservation is a consequence of the constraint geometry, not an independent postulate.

6.2. Pass-1 Exactness and the Role of Height

Theorem 1 establishes that the speed field requires only $z\left(\lambda \right)$ , regardless of parametrization. For all eight families in Table 3, $z\left(\lambda \right)$ is a simple closed-form expression, making Pass 1 always exactly solvable. For the cubic track, this means the full rail loading map $N\left(s\right)$ —including the inward/outward loading transition and the exact $N=0$ crossing—is available without any numerical integration, even as the curvature changes sign continuously through the inflection point.

6.3. Continuous Frame through Inflection Points

A key feature of GQM, demonstrated in Example 2, is the use of an algebraically signed curvature $\kappa \left(s\right)$ and a globally oriented normal $\widehat{n}\left(s\right)$ . In standard Frenet-frame implementations, the normal is constrained to point toward the center of curvature, forcing a 180˚ flip when the curve passes an inflection point. This flip produces an artificial jump in the computed normal force $N\left(s\right)$ that has no physical basis. By defining $\widehat{n}$ via a fixed counter-clockwise rotation of the tangent (Equation (7)), GQM avoids this discontinuity entirely: $\kappa \left(s\right)$ passes smoothly through zero at the inflection, $\widehat{n}\left(s\right)$ remains continuous, and $N\left(s\right)$ transitions continuously between inward and outward loading without any special treatment at the inflection boundary.

6.4. Inescapability of Non-Elementary Functions in Pass 2

By Liouville’s theorem on integration in finite terms, as formalized by Ritt [13], the time integral ${\displaystyle \int \text{d}s/\sqrt{c-2gz\left(s\right)}}$ is non-elementary whenever $z\left(s\right)$ introduces algebraic functions of degree exceeding two. (The non-elementary character of elliptic integrals is the classical special case, with the general non-elementarity criterion formalized by Liouville and Ritt extending it to the broader setting used here.) For circular arcs this is $K\left(\text{sin}\left({\theta }_0/2\right)\right)$ ; for the cubic $z=x^3$ the height function introduces a quartic algebraic relationship through the arc-length element. No reformulation can circumvent this: GQM confines the non-elementary computation to a single scalar quadrature in Pass 2, keeping all spatial quantities (Pass 1) exactly elementary.

6.5. The Pendulum as Validator, Not Primary Target

For the pendulum, the motion is governed by a well-known equation whose long-time behavior can be tracked by high-order integrators to high but not unlimited accuracy. GQM reconstructs the motion using exact half-period tiling based on high-accuracy spatial quadrature, reducing phase and energy drift by many orders of magnitude compared with direct ODE integration (Table 4). Both RK4 and Störmer-Verlet accumulate phase error that grows as $O\left(h^p T_{\text{final}}\right)$ ; as shown in (32), GQM’s total pipeline phase error after $N$ periods is bounded, in time-domain terms, by $N{\epsilon }_{\tau }+{\epsilon }_{\text{interp}}\left(\delta t\right)$ , where ${\epsilon }_{\tau }\approx {10}^{-13}\text{-}{10}^{-15}$ and the interpolation term does not grow with $N$ [15]; the corresponding angular-equivalent bound, obtained by multiplying through by the local angular rate $\dot{\theta }$ , is correspondingly smaller than RK4’s angular drift by many orders of magnitude. The pendulum thus serves as a validation case, not the primary target of the method.

For any other smooth planar curve—ellipses, spirals, catenaries, cycloids, arbitrary surface profiles—no closed-form trajectory formula exists or can exist (Liouville’s theorem [13]). In this universally applicable regime, among the methods compared in Table 7—direct ODE integration, symplectic integration, perturbation/series expansion, and special-function solutions—GQM is, to the authors’ knowledge, the only one that simultaneously provides exact Pass-1 spatial quantities, machine-precision time-of-flight by a single quadrature, and phase-error bounded at $N{\epsilon }_{\tau }$ at arbitrarily large time. This claim is scoped to the comparison performed here and to the class of conservative planar-curve problems satisfying assumptions (A1) - (A4) of Section 2.3; it is not asserted as an exhaustive survey of the broader numerical-methods literature. Table 7 summarises.

Table 7. Comparison of solution methods. “Exact” = machine precision in one evaluation. “N-exact, drift $\le N{\epsilon }_{\tau }$ ” = machine precision per period; accumulated phase error bounded by $N{\epsilon }_{\tau }+{\epsilon }_{\text{interp}}\left(\delta t\right)$ , with ${\epsilon }_{\tau }\approx {10}^{-14}$ and ${\epsilon }_{\text{interp}}$ non-accumulating across periods (see (32)). “Drifts” = phase error grows as $O\left(h^p{T}_{\text{final}}\right)$ . The bold row is the present method. ^†marks structurally empty columns.

Method	Circle/pendulum only		Arbitrary smooth curve
	$V^2\left(s\right)$ , $N\left(s\right)$	$\theta \left(t\right)$ , large $t$	$V^2\left(s\right)$ , $N\left(s\right)$	$\lambda \left(t\right)$ , large $t$
Jacobi/special functions	Exact	Exact, no drift	†None exists	†None exists
Direct ODE (RK4)	N-exact	N-exact, drifts	N-exact	N-exact, drifts
Symplectic (Störmer-Verlet)	N-exact	Energy bounded; phase drifts	N-exact	Energy bounded; phase drifts
Perturbation/series	Approx.	Small amplitude only	Approx.	Restricted validity
Spatial method (GQM)	Exact (P1)	N-exact, drift $\le N{\epsilon }_{\tau }$	Exact (P1)	N-exact, drift $\le N{\epsilon }_{\tau }$

6.6. Phase Stability at Large Time

The analysis in this section applies to conservative (periodic) systems, for which Pass 3 exact half-period tiling is applicable. For dissipative systems (e.g. Coulomb friction, Example 2), the motion is non-periodic and Pass 3 reduces to direct interpolation on the Pass-2 table without tiling; the phase-stability comparison below does not apply to that case.

A $p$ -th order explicit integrator accumulates a local truncation error of $O\left(h^{p+1}\right)$ per step, giving a global phase error of $O\left(h^p\right)$ per period, growing to

${\epsilon }_{\text{phase}}\left(T_{\text{final}}\right)=O\left(h^p{T}_{\text{final}}\right)$ (31)

over the full integration. Symplectic integrators conserve a shadow Hamiltonian exactly, suppressing secular energy growth, but the shadow Hamiltonian differs from the true one by $O\left(h^p\right)$ , giving a period error of $O\left(h^p\right)$ per cycle and accumulated phase error of $O\left(h^p{T}_{\text{final}}\right)$ —identical in scaling to RK4 [15].

GQM bounds drift via two contributions that must be accounted for separately. The first is the quadrature error: the half-period $\tau$ is computed once by adaptive Gauss-Kronrod quadrature [14] to absolute precision ${\epsilon }_{\tau }\approx {10}^{-13}\text{-}{10}^{-15}$ , and exact half-period tiling then contributes a phase error of at most $N{\epsilon }_{\tau }$ after $N$ periods. The second is the Pass-3 inversion error: $s\left(t\right)$ is recovered by interpolation on the precomputed Pass-2 table $\left\{t_i,s_i\right\}$ , and the interpolation error ${\epsilon }_{\text{interp}}$ depends on the table resolution $\delta t$ and the chosen scheme. For a cubic-spline interpolant on a uniform table with spacing $\delta t$ , the inversion error is $O\left(\delta t^4\right)$ per query, independent of $N$ . The total pipeline phase error at $N$ periods is therefore bounded by

${\epsilon }_{\text{phase}}^{\text{GQM}}\le N{\epsilon }_{\tau }+{\epsilon }_{\text{interp}}\left(\delta t\right),$ (32)

where the second term does not grow with $N$ . In the numerical experiments of Table 4, a cubic-spline table with $\delta t={10}^{-4}$ s was used. For the $L=1\text{m}$ pendulum ($T_{\text{max}}=1000\text{s}$ , period $\tau \approx 2.20\text{s}$ , $N\approx 454$ periods), this gives ${\epsilon }_{\text{interp}}\approx {10}^{-14}$ and a tiling contribution $N{\epsilon }_{\tau }\approx 454\times {10}^{-15}\approx 5\times {10}^{-13}$ —many orders of magnitude below both RK4 and symplectic integrators, with the interpolation error non-accumulating across periods. To illustrate the scaling advantage at larger horizons: for a hypothetical run of $N={10}^6$ periods with ${\epsilon }_{\tau }\approx {10}^{-15}$ , the bound gives $N{\epsilon }_{\tau }\approx {10}^{-9}$ , still far below ODE-integrator drift of $O\left(h^p{T}_{\text{final}}\right)$ . This bound applies to the class of conservative planar-curve problems addressed by this paper (assumptions (A1) - (A4)); within this class, and among the standard time-stepping and series-based alternatives compared in Table 7, no competing approach is known to the authors that achieves phase stability of this form for general (non-circular) planar curves.

6.7. Limitations

The present formulation assumes:

1) A planar (2D) curve;

2) A particle (point mass), though the rolling extension admits rigid bodies with a single rolling degree of freedom;

3) Gravity as the primary potential force (the rotating-frame kernel extends this to centrifugal potential, and any conservative force can be accommodated by replacing the gravitational kernel $\mathcal{K}$ as described in Section 3).

Extension to space curves (3D) requires the full Darboux frame and is deferred to future work.

Beyond the dimensional and force-model restrictions above, the formulation also relies on the regularity assumptions (A1) - (A2) of Section 2.3, which may fail or require modification in the following cases:

Non-smooth tracks and cusps.

At a cusp or any point where $\text{d}r/\text{d}\lambda \to 0$ or the tangent direction is discontinuous, the curve fails to be $C^1$ , the unit-speed parametrization (4) cannot be extended across the point, and the curvature (8) is undefined there. GQM as formulated requires a piecewise treatment in this case: each smooth segment is parametrized and solved independently by the three-pass pipeline, and the segments are joined by a kinematic hand-off analogous to the discrete stitching of Section 2.1 (continuity of position and speed, but not necessarily of tangent direction, across the cusp). This is a natural but currently unimplemented extension of the piecewise-linear limiting case already used to motivate the method.

Self-intersecting curves.

The arc-length map $s\mapsto \left(x\left(s\right),z\left(s\right)\right)$ is one-to-one by construction, but if the underlying curve self-intersects, several distinct values of $s$ correspond to the same point $\left(x,z\right)$ in the plane. GQM’s spatial quantities (Pass 1) remain valid along the single-valued arc-length path, since all computations are functions of $s$ rather than of $\left(x,z\right)$ directly. However, any query that requires identifying where on the physical curve the particle currently is—e.g. for collision detection with another branch of the same curve—requires additional bookkeeping beyond the arc-length state, since $\left(x,z\right)$ alone no longer determines $s$ uniquely. The present formulation implicitly assumes a simple (non-self-intersecting) curve over the domain of integration; self-intersections are not treated.

Separatrix-like and asymptotic motions.

Pass 3 tiling (18) assumes a finite half-period $\tau$ . For a trajectory approaching an unstable equilibrium asymptotically (e.g. a pendulum released exactly at the inverted position, or more generally any orbit asymptotic to a separatrix), $\tau \to \infty$ and the Pass-2 quadrature (17) converges arbitrarily slowly near the equilibrium, since the integrand’s singularity there is no longer the generic ${\left(s_{\star }-s\right)}^{-1/2}$ algebraic type assumed in Section 2.7 but a logarithmically divergent one. In this boundary regime the finite-quadrature-cost guarantee of Pass 2 and the bounded-tiling guarantee of Pass 3 (Section 6.6) no longer hold without separate asymptotic treatment (e.g. analytic matching near the equilibrium), which is not developed in the present work.

Loss of monotonicity of t(s).

Pass 3 explicitly assumes (Section 2.7) that $t\left(s\right)$ is strictly monotone on each arc between turning points, which holds whenever $V\left(s\right)$ maintains a single sign on that arc. This can fail in two ways: 1) for dissipative kernels with multiple sign changes of $N$ (Section 3.2), if the restart procedure of Section 3.2 is not applied correctly the computed $V^2\left(s\right)$ may become spuriously negative or oscillate near a crossing, corrupting monotonicity of $t\left(s\right)$ ; and 2) more fundamentally, if the force model itself permits $V$ to reverse sign without $V=0$ being reached exactly (not possible for the conservative and friction/drag kernels of Table 2, but conceivable for a more general velocity-dependent kernel $\mathcal{A}\left(s;V\right)$ not considered here). In either case the Pass-3 interpolation/inversion step becomes ill-posed and the GQM pipeline as stated does not apply without modification.

7. Conclusions

The Geometric Quadrature Method (GQM) reformulates constrained particle dynamics on arbitrary planar curves as a three-pass spatial pipeline, delivering capabilities unavailable from standard ODE-based methods.

Pass 1 yields the speed field and normal contact force in exact closed form for any curve with an analytic height function—no quadrature, no approximation. This is the core structural advantage: the entire rail loading map, including force sign transitions and liftoff conditions, is accessible algebraically. Pass 2 confines the irreducibly non-elementary time-of-flight computation to a single adaptive quadrature, evaluated to machine precision. Pass 3 recovers phase-stable long-time trajectories for conservative (periodic) systems by exact half-period tiling, bounding phase error at $N{\epsilon }_{\tau }$ independently of $T_{\text{final}}$ ; for dissipative kernels the motion is non-periodic and Pass 3 reduces to direct table inversion without tiling.

A central element enabling the method’s extensibility is the speed-field kernel $\mathcal{A}\left(s\right)$ , which encodes the physical force model in a modular, plug-and-play form. The kernel specifies how each arc-length element contributes to the total speed change—in exact analogy with the convolution kernel of the discrete piecewise-linear limit (Section 2.1). Six kernels are derived covering frictionless sliding, rolling with inertia, Coulomb friction, quadratic drag, surface tension, and rotating frames. A complementary dimensional reduction table (Table 3) reduces eight standard curve families—Cartesian graphs, circles, ellipses, parabolas, logarithmic spirals, catenaries, cycloids, and implicit curves—to this single framework by providing the arc-length element and curvature for each geometry. Importantly, any conservative force field can be accommodated: replacing the gravitational term in $\mathcal{A}$ with $-\nabla U\cdot \widehat{t}$ for an arbitrary potential $U\left(x,z\right)$ yields the generalized Pass-1 result $V^2\left(s\right)=V_0^2+2\left[U_0-U\left(s\right)\right]$ , with Pass 2 and Pass 3 proceeding unchanged.

The method’s position is clear. For the simple pendulum, the phase-stability benchmark shows that GQM reduces accumulated phase and energy drift by many orders of magnitude compared with all standard time-stepping schemes—explicit and symplectic alike—over long simulations (see Equation (32) and Table 4). For every other smooth planar curve—ellipses, spirals, catenaries, cycloids, cubic profiles, and arbitrary numerically defined tracks—Liouville’s theorem [13] guarantees that no closed-form trajectory exists. In this vast, practically relevant regime, and among the methods surveyed in Table 7, GQM is the only approach we are aware of that provides both exact spatial quantities and phase-stable trajectories, within the conservative, planar, point-mass setting of assumptions (A1) - (A4).

The cubic track example further demonstrates a geometric robustness advantage unique to GQM: the use of an algebraically signed normal and a continuously signed curvature eliminates the frame-flip artefacts that conventional Frenet-frame solvers produce at inflection points. This makes GQM directly applicable to realistic engineering tracks that combine convex and concave segments.

Practical applications include rail loading prediction, long-time phase-stable simulation of non-circular oscillators, parametric sweeps over curve families in design optimization, and any setting where exact spatial force maps are needed without full trajectory integration. The kernel substitution principle of Section 3 opens further domains: replacing the gravitational term with a central-force gradient $-\left(\mu /r^2\right)\text{d}r/\text{d}s$ recovers the vis-viva equation in Pass 1 and the Kepler time-of-flight integral in Pass 2, extending naturally to perturbed non-Keplerian orbits. For geometries that evolve slowly relative to the particle’s oscillation frequency—vibrating rails, flexible beams—the kernel can be updated quasi-statically each period, preserving Pass-1 exactness in the adiabatic regime.

Future work will extend the framework to 3D space curves via the Darboux frame, and to impact and rebound modeling.

Data Availability

The Python script GQM_Plots.py used to generate all numerical results and figures is provided as supplementary material. No other datasets were generated or analysed during the current study.

Appendix

Curvature Formulas for Standard Parametrizations

Cartesian graph $z=\phi \left(x\right)$ :

$\kappa =\frac{{\phi }^{″}}{{\left(1+{{\phi }^{\prime }}^2\right)}^{3/2}}.$ (33)

Parametric curve $r\left(\lambda \right)=\left(x\left(\lambda \right),z\left(\lambda \right)\right)$ :

$\kappa =\frac{\dot{x}\ddot{z}-\dot{z}\ddot{x}}{{\left({\dot{x}}^2+{\dot{z}}^2\right)}^{3/2}},$ (34)

where dots denote $\text{d}/\text{d}\lambda$ .

Polar curve $r=f\left(\theta \right)$ :

$\kappa =\frac{f^2+2{f^{\prime }}^2-f{f}^{″}}{{\left(f^2+{f^{\prime }}^2\right)}^{3/2}}.$ (35)

Implicit curve $F\left(x,z\right)=0$ :

$\kappa =\frac{F_{xx}{F}_z^2-2F_{xz}{F}_x{F}_z+F_{zz}{F}_x^2}{{\left|\nabla F\right|}^3}.$ (36)