Universal Polynomial Invariants of Nonholonomic Leapfrog

1. Introduction

Leapfrog is a numerical method for approximating solutions of mechanical differential equations through staggered position and velocity updates. The nonholonomic version considered here, which also enforces a velocity constraint, was proposed by McLachlan and Perlmutter in 2006 [1] and revisited through a geometric construction by Ferraro, Iglesias, and Martín de Diego in 2008 [2]; later extensions and convergence results appeared in 2015 [3]. In 2020, Modin and Verdier [4] reported favorable leapfrog energy behavior in a perturbed continuously variable transmission (CVT) experiment that their reversibility argument did not explain. In 2023, Anahory Simoes et al. [5] revisited that model and observed favorable energy behavior for a Newmark composition method. These numerical observations motivate our study of the conservation laws built into the reflected leapfrog update.

1.1. The CVT Model and Its Energy Observation

The continuously variable transmission (CVT) model couples a spring subsystem, called the passenger, to a driver that changes its allowed velocity direction. With normalized masses and spring stiffnesses, passenger position $x=(x_1,x_2)$ and driver coordinate $\xi$ satisfy (Equation (12) [4])

$$\begin{array}{cccc}\hfill {\ddot{x}}_1& =-x_1+\lambda ,\hfill & \hfill {\ddot{x}}_2& =-x_2+\lambda sin\xi ,\hfill \\ \hfill \ddot{\xi }& =-V^{\prime }\left(\xi \right),\hfill & \hfill {\dot{x}}_1+sin\xi {\dot{x}}_2& =0.\hfill \end{array}$$

(1)

Dots denote time derivatives. The potential V determines the independent driver motion; the constraint force $\lambda (1,sin\xi )$ enforces the allowed passenger velocity direction. With $p=\dot{x}$, the passenger energy is

$$E=\frac{1}{2}{(\parallel x\parallel }^2+{\parallel p\parallel }^2),$$

the sum of kinetic and spring potential energies. The constraint force does no work because $\lambda (p_1+sin\xi p_2)=0$, so the continuous model conserves E. The numerical question is how well leapfrog retains this conservation. For the reconstruction below we use $V\left(\xi \right)=1-cos\xi -\frac{1}{4}sin\left(2\xi \right)$; the caption records its distinction from the printed setup.

Modin and Verdier (Figure 6 [4]) compared numerical methods on this model. In the rotating, perturbed case, leapfrog maintained small errors in the energy of the spring subsystem over the reported interval $0\le t\le 3000$ with $h=0.1$. Their reversibility argument did not explain this behavior, and they posed its explanation as an open problem (Section 6.1(2) [4]). Figure 1 shows our independent reconstruction of the qualitative contrast: leapfrog’s passenger-energy error stays small over the simulated interval while the midpoint discrete Lagrange–d’Alembert (DLA) error grows.

Anahory Simoes et al. (Section 5.2 [5]) later revisited the perturbed CVT and reported similarly favorable energy behavior for a Newmark composition method, alongside comparisons with leapfrog. Its explanation also remained open in that study. We analyze the reflected leapfrog update below; extending the classification to the Newmark composition would require a separate analysis.

Figure 2 summarizes the methods and numerical observations motivating our study.

1.2. From the Observation to a Classification Question

Figure 3 follows our investigation from the authors’ question to the estimate that remains open. Modin and Verdier ask what explains leapfrog’s favorable passenger-energy behavior in the rotating, perturbed CVT experiment (Section 6.1(2) [4]). Earlier reflection methods and convergence results provide the starting framework [1,2,3], but Modin and Verdier’s reversibility argument does not cover this case. These are the starting question and prior foundations in steps 1–2.

Our classification (step 3) isolates what the update conserves independently of the driver. The CVT velocity constraint in (1) has projector

$$P\left(\xi \right)=a\left(\xi \right)a{\left(\xi \right)}^T,a\left(\xi \right)={\displaystyle \frac{{(-sin\xi ,1)}^T}{\sqrt{1+{sin}^2\xi }}}.$$

We retain the reflected spring update and allow P to range over all planar rank-one orthogonal projectors, a larger family than the directions traced by this CVT driver. For planar rank-one projectors and $0<h<2$, every common polynomial invariant is a polynomial in the modified quadratic energy $J_h$, except at $h=\sqrt{2}$, where an additional generator $W_h^2$ appears. Special oscillator step sizes already occur in [1]; our contribution is completeness for the entire update family.

The energy identity (step 4) connects this classification to the original question. Writing $C_h=J_h-E$ for the correction in (19), conservation of $J_h$ gives $E_n-E_0=C_0-C_n$, where $C_n$ is the correction along the numerical trajectory. At $h=0.1$, no independent common polynomial supplies additional information beyond $J_h$. Our next target (step 5) is therefore to use the actual CVT driver to bound the variation of $C_n$ over $0\le nh\le 3000$. The identity is established here; the small bound needed to explain the observation is not. Section 2, Section 3 and Section 4 develop the classification in step 3 of Figure 3; Section 5 establishes the energy identity in step 4 and states the remaining target in step 5.

Next, in Section 2, we write down how the method advances the spring and identify a quantity that remains unchanged after each step.

2. The Update and the Classification Problem

To begin proving the classification in step 3 of Figure 3, we specify the update family, identify its conserved quadratic energy, and state the full result in Theorem 1. We begin in Section 2.1 by writing down one step of the numerical method: how it calculates the next position and velocity when the allowed direction of motion can change.

2.1. The Reflected Spring Update

Let $x,u\in {\mathbb{R}}^2$, where x is position and u is the incoming half-step velocity stored by the method. The stored velocity need not lie in the currently allowed direction. The spring stiffness and mass are normalized to one. For an orthogonal rank-one projector P, define $F_{P,h}(x,u)=(x^{\prime },u^{\prime })$ by

$$u^{\prime }=(2P-\mathrm{I})u-hPx,x^{\prime }=x+h{u}^{\prime }.$$

(2)

The projector selects the allowed velocity direction. The map $2P-\mathrm{I}$ reflects across that direction. The synchronized velocity is $p=(u+u^{\prime })/2=P(u-hx/2)$, which satisfies the current constraint. We fix $0<h<2$ and allow a different P at every step.

With the update specified, Section 2.2 asks what stays unchanged as we repeat it. We first find an energy-like quantity that keeps the same value, even when the allowed direction changes, and then state which other conserved polynomials we aim to identify.

2.2. The Modified Energy and the Classification Question

We will use the quadratic quantities

$$J_h={\displaystyle \frac{{\parallel x\parallel }^2+{\parallel u\parallel }^2-h\langle x,u\rangle }{2}},W_h=h\langle x,u\rangle -{\parallel u\parallel }^2.$$

(3)

The first is conserved for every projector; the second will identify the exceptional step. The proof below uses the norm-preserving reflection underlying [2]. We include the calculation for this spring update to make its arbitrary-projector scope explicit.

Proposition 1

(Modified energy). The update (2) preserves $J_h$ for every orthogonal projector P and every real h. Consequently $J_h$ remains constant along arbitrary projector sequences with fixed h.

Proof. 

Write $v=u-hx/2$ and $R_P=2P-\mathrm{I}$. Then $u^{\prime }+hx/2=R_{P}v$. Since $R_P$ is orthogonal,

$${\parallel u^{\prime }\parallel }^2-{\parallel u\parallel }^2=-h\langle x,u^{\prime }+u\rangle .$$

(4)

Expanding $J_h(x+h{u}^{\prime },u^{\prime })-J_h(x,u)$ and applying (4) gives zero at each step. □

We now ask whether any independent polynomial is conserved under the same freedom to change P.

Definition 1

(Common polynomial invariant). Fix $0<h<2$. A polynomial $f(x,u)$ is a common invariant if

$$f\left(F_{P,h}(x,u)\right)=f(x,u)$$

for every $x,u\in {\mathbb{R}}^2$ and every rank-one orthogonal projector P. We denote the algebra of these polynomials by ${\mathcal{A}}_h$.

Theorem 1

(Complete classification). For the family (2) with $x,u\in {\mathbb{R}}^2$ and $0<h<2$,

$${\mathcal{A}}_h=\left\{\begin{array}{cc}\mathbb{R}\left[J_h\right],\hfill & h\ne \sqrt{2},\hfill \\ \mathbb{R}[J_h,W_h^2],\hfill & h=\sqrt{2}.\hfill \end{array}\right.$$

(5)

At $h=\sqrt{2}$, the generators $J_h$ and $W_h^2$ are algebraically independent.

In other words, every common polynomial invariant is a polynomial in the displayed generators, with no restriction on degree. At the exceptional step, the additional generator has degree four; the space of homogeneous quadratic invariants is still one-dimensional. Next, in Section 3, we show where the extra quartic comes from: at $h=\sqrt{2}$, the update exchanges two squared lengths, producing another conserved quantity.

3. Finding the Exceptional Generator

The modified energy supplies one polynomial generator. We now construct the additional quartic at $h=\sqrt{2}$ by changing coordinates; completeness will require a separate argument.

Quartic integrals already occur in discrete nonholonomic mechanics: Fedorov and Zenkov (Theorem 6.2 [6]) give an integral quartic in quaternion variables for the discrete Suslov system. Here the question is conservation under every projector choice and the completeness of the resulting algebra. The extra generator comes from an exchange of lengths.

3.1. A Change of Variables Exchanges Lengths

Introduce

$$y=hu-x.$$

(6)

Substitution into (2) gives, for every nonzero h,

$$x^{\prime }=(2-h^2)Px+R_{P}y,y^{\prime }=-x.$$

(7)

At $h=\sqrt{2}$ the first term vanishes. The remaining reflection preserves the length of y, while the sign change preserves the length of x. Therefore

$${\parallel x^{\prime }\parallel }^2={\parallel y\parallel }^2,{\parallel y^{\prime }\parallel }^2={\parallel x\parallel }^2.$$

(8)

The squared lengths exchange values. Their sum and product remain unchanged, as illustrated in Figure 4.

3.2. The Quartic Invariant

The length exchange yields an exact conservation law, even beyond the planar setting of our classification.

Proposition 2

(Resonant conservation in any dimension). Let $x,u$ belong to a real inner-product space, let $h^2=2$, and use (2) with any orthogonal projector P. Then

$$W_h(x^{\prime },u^{\prime })=-W_h(x,u).$$

Consequently $W_h^2$ and ${\parallel x\parallel }^2{\parallel hu-x\parallel }^2$ are conserved along every sequence of such updates, with fixed h.

Proof. 

The calculation (7) and the length exchange (8) use only $P^2=P$ and $P^{*}=P$, so they hold in any dimension. At $h^2=2$,

$${\parallel x\parallel }^2+{\parallel y\parallel }^2=4J_h,{\parallel x\parallel }^2-{\parallel y\parallel }^2=2W_h.$$

(9)

The difference changes sign under the exchange, and the product is unchanged. In particular, ${\parallel x\parallel }^2{\parallel y\parallel }^2=4J_h^2-W_h^2$. □

For example, take a unit vector e. The states $(x,u)=(e,0)$ and $(0,e)$ both have $J_h=1/2$, whereas $W_h^2$ equals 0 and 1, respectively. The quartic therefore distinguishes states that the quadratic energy does not. Next, in Section 4, we show that these quantities account for every common polynomial invariant. We do this by studying the transformations built from repeated updates and what those transformations leave unchanged.

4. Proving that the Generators Are Complete

Finding conserved quantities gives only part of the classification. We now show that every common polynomial invariant belongs to the algebra generated by them. First we identify the group generated by the updates; then we determine its invariant polynomials.

4.1. Identifying the Group Generated by the Updates

A common invariant must survive every composition of the updates. The group generated by all constraint directions therefore determines what an invariant can distinguish. We work in the coordinates $z={(x_1,x_2,u_1,u_2)}^T$ and write

$$A_h=\left(\begin{array}{cccc}1-h^2& 0& h& 0\\ 0& 1& 0& -h\\ -h& 0& 1& 0\\ 0& 0& 0& -1\end{array}\right),G_h=\left(\begin{array}{cc}{I}_2& -h{I}_2/2\\ -h{I}_2/2& I_2\end{array}\right),Z=\left(\begin{array}{cccc}0& 1& 0& 0\\ -1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& -1& 0\end{array}\right).$$

(10)

Here $A_h$ is the update for $P=diag(1,0)$, and $D\left(\alpha \right)=e^{\alpha Z}$ rotates position and velocity together. Every rank-one orthogonal projector therefore gives an update $D\left(\alpha \right)A_{h}D(-\alpha )$. Direct multiplication gives

$$A_h^T{G}_h{A}_h=G_h,Z^T{G}_h+G_{h}Z=0,det{A}_h=-1,z^T{G}_{h}z=2J_h\left(z\right).$$

(11)

The eigenvalues of $G_h$ are $1\pm h/2$, each twice. Thus $G_h$ is positive definite on the interval under consideration.

Let $H_h$ be the closed group generated by the projector updates. For $S_h=G_h^{1/2}$, the group $S_h{H}_h{S}_h^{-1}$ is a closed subgroup of $O\left(4\right)$. In these coordinates, $2J_h$ becomes ordinary squared length. Here $O\left(4\right)$ denotes all length-preserving linear maps in four dimensions, and $SO\left(4\right)$ its determinant-one subgroup. The notation ${\mathbb{T}}^2$ describes two independent plane rotations; $C_4$ is the cyclic group of order four. At resonance their semidirect product records how an order-four update rearranges those rotations.

Theorem 2.

For $0<h<2$,

$$S_h{H}_h{S}_h^{-1}=O\left(4\right)\mathrm{if}{h}^2\ne 2.$$

If $h^2=2$, the group is isomorphic to ${\mathbb{T}}^2⋊C_4$. Its identity component independently rotates two orthogonal planes, and a generator of $C_4$ exchanges those planes.

Proof. Step 1: Recover simultaneous rotations. We first show that $H_h$ contains all simultaneous rotations $D\left(\alpha \right)$. In this proof write $A=A_h$, $G=G_h$, and

$$X_j=A^{j}Z{A}^{-j},Y_j=X_j-X_{j+1},d=h^2(h^2-4).$$

The curve $D\left(\alpha \right)AD(-\alpha )A^{-1}$ lies in $H_h$ and starts at the identity. Its derivative is $Y_0$, so $Y_0$ belongs to the Lie algebra of $H_h$. Conjugating by A puts every $Y_j$ in that Lie algebra. The matrices in (10) satisfy

$$dZ=(d+1)(Y_0+Y_1)-(Y_2+Y_3).$$

(12)

Since $d\ne 0$, this puts Z in the Lie algebra. Consequently

$$H_h=\overline{\langle A,D\left(\alpha \right):\alpha \in \mathbb{R}\rangle }.$$

(13)

Step 2: Generate the group away from resonance. Suppose first that $h^2\ne 2$. Here we seek matrix generators for the Lie algebra, rather than generators of conserved polynomials. The Lie algebra contains the six matrices

$$X_0,X_1,X_2,X_3,[X_2,X_0],[X_3,X_1],[B,C]=BC-CB.$$

(14)

For each matrix, retain its entries in positions $(1,1),(1,2),(1,4),(2,1),(2,2),(3,2)$, in that order, and form the six-column matrix using (14). Its determinant is

$$2h^8{(h^2-2)}^4.$$

(15)

Equations (12) and (15) are finite matrix identities; (10) and the specified entries give all data needed to check them. The six matrices are independent and G-skew, so they span the six-dimensional Lie algebra $\{B:B^{T}G+GB=0\}$. After conjugation by $S_h$, the identity component is therefore $SO\left(4\right)$. Since $detA=-1$, the full group is $O\left(4\right)$.

Step 3: Identify the exceptional group. Now suppose $h^2=2$, and put $K=AZ{A}^{-1}$. Direct multiplication gives

$$Z^2=K^2=-I_4,ZK=KZ,AK{A}^{-1}=-Z,A^4=I_4.$$

(16)

The matrix

$$T=ZK=\left(\begin{array}{cc}{I}_2& 0\\ h{I}_2& -I_2\end{array}\right)$$

(17)

is a G-self-adjoint involution. Its eigenspaces $V_{+}$ and $V_{-}$ have dimension two and are G-orthogonal. On $V_{+}$, $K=-Z$; on $V_{-}$, $K=Z$. Hence the subgroup

$$\mathcal{T}=\{e^{sZ+tK}:s,t\in \mathbb{R}\}$$

independently rotates these planes and is a two-dimensional torus. Conjugation by A maps $(Z,K)$ to $(K,-Z)$, so it normalizes $\mathcal{T}$. No element $A^j$, $j=1,2,3$, lies in $\mathcal{T}$: its conjugation sends Z to $K,-Z,-K$, respectively, whereas all elements of $\mathcal{T}$ commute with Z. Thus $\mathcal{T}\langle A\rangle$ is a compact group isomorphic to ${\mathbb{T}}^2⋊C_4$. By (13) it is exactly $H_h$. Finally, $AT{A}^{-1}=-T$, so A exchanges $V_{+}$ and $V_{-}$. Conjugation by $S_h$ makes these planes Euclidean-orthogonal. □

4.2. Recovering All Invariant Polynomials

The group description now completes the classification: away from resonance, only the total squared length matters; at resonance, the two plane lengths supply the additional generator.

Proof 

(Proof of Theorem 1). A polynomial preserved by each invertible update is also preserved by its inverse, by every finite product, and by the closure of those products. Thus the common polynomial invariants are precisely the $H_h$-invariant polynomials.

If $h^2\ne 2$, Theorem 2 reduces the question to $O\left(4\right)$. An invariant polynomial is constant on every sphere. Its restriction to a coordinate axis is an even polynomial, say $q\left(t^2\right)$; invariance then gives $q(\parallel S_{h}z{\parallel }^2)=q\left(2J_h\left(z\right)\right)$ everywhere. This proves that the invariant algebra is $\mathbb{R}\left[J_h\right]$.

At $h^2=2$, let $r_{+}$ and $r_{-}$ be the squared $G_h$-lengths of the components in the two planes from (17):

$$r_{\pm }=\frac{1}{2}{z}^T{G}_h(I_4\pm T)z,r_{+}+r_{-}=2J_h,r_{+}-r_{-}=W_h.$$

(18)

The last equality follows from $G_{h}T=\left(\begin{array}{cc}0& h{I}_2/2\\ h{I}_2/2& -I_2\end{array}\right)$. In the reflection coordinates $y=hu-x$, these are simply $r_{+}={\parallel x\parallel }^2/2$ and $r_{-}={\parallel y\parallel }^2/2$. Independent rotations on the two planes force every invariant polynomial to belong to $\mathbb{R}[r_{+},r_{-}]$. Indeed, in complex coordinates on each plane, rotation invariance forces each monomial to contain equal powers of the coordinate and its conjugate. The generator $A_h$ exchanges the two radii. The full invariant algebra is therefore the algebra of symmetric polynomials in them:

$$\mathbb{R}[r_{+}+r_{-},r_{+}{r}_{-}]=\mathbb{R}[J_h,W_h^2],r_{+}{r}_{-}=J_h^2-\frac{1}{4}{W}_h^2.$$

There is no polynomial relation between these two generators. The radii range independently over the nonnegative quadrant, and the map $(r_{+},r_{-})\mapsto (J_h,W_h^2)$ has nonzero Jacobian when $r_{+}\ne r_{-}$. Its image therefore contains an open set, on which no nonzero polynomial can vanish. This completes the classification. □

Next, in Section 5, we relate these conserved quantities to the spring subsystem’s physical energy and identify the estimate still needed to explain the observed behavior at $h=0.1$.

5. From Modified Energy to the Remaining CVT Estimate

The common polynomial classification is now complete. Steps 4–5 of Figure 3 concern a different task: connecting that conservation to physical energy and identifying the estimate needed for the observed CVT behavior.

5.1. The Physical-Energy Identity

To relate the conserved quantity to physical energy, define the synchronized velocity $p=P(u-hx/2)$ and the remaining numerical component $r=(\mathrm{I}-P)(u-hx/2)$. Orthogonality gives

$$J_h=E+{\displaystyle \frac{{\parallel r\parallel }^2}{2}}-{\displaystyle \frac{h^2{\parallel x\parallel }^2}{8}},E={\displaystyle \frac{{\parallel x\parallel }^2+{\parallel p\parallel }^2}{2}}.$$

(19)

Thus conservation of $J_h$ does not by itself imply conservation of physical energy. The energy-error estimate depends on the variation of the combined correction. The additional invariant at $h=\sqrt{2}$ is an exact property of a particular step size, not a small-step accuracy result.

5.2. Future Work: The CVT Estimate at $h=0.1$

At step n, use the projector $P_n$ at position $x_n$ and define

$$\begin{array}{cccc}\hfill p_n& =P_n(u_n-h{x}_n/2),\hfill & \hfill r_n& =(\mathrm{I}-P_n)(u_n-h{x}_n/2),\hfill \\ \hfill E_n& =\frac{1}{2}(\parallel x_n{\parallel }^2+\parallel p_n{\parallel }^2),\hfill & \hfill C_n& =\frac{1}{2}\parallel r_n{\parallel }^2-\frac{h^2}{8}{\parallel x_n\parallel }^2.\hfill \end{array}$$

Here $p_n$ is the synchronized passenger velocity; $u_n$ remains the stored half-step velocity. Equation (19) and conservation of $J_h$ give $E_n-E_0=C_0-C_n$. Our proposed target for the reported experiment is a small quantitative bound on

$$\underset{0\le nh\le 3000}{max}|C_n-C_0|,h=0.1,$$

using the actual CVT driver and initial data. This estimate is not proved here. It concerns the variation of the combined correction; neither term must vanish separately. An estimate uniform over all n would establish the stronger all-time result.

The next step is to use the ordered projector sequence $P_n=P\left({\xi }_n\right)$ generated by the numerical CVT driver to estimate $C_n-C_0$. One route is a direct bound that accounts for cancellation between steps. Another is to seek a conserved or nearly conserved quantity depending on the driver’s phase as well as on $(x,u)$, and quantify its relation to $C_n$. Theorem 1 does not exclude such driver-dependent quantities. Either route must produce a bound small enough to account for the observed energy error, rather than merely show that the trajectory stays bounded.

6. Conclusions

We have answered which polynomial conservation laws are shared by every planar rank-one projector update of nonholonomic leapfrog. For fixed $0<h<2$, the modified energy $J_h$ generates all common polynomial invariants, except at $h=\sqrt{2}$, where the additional generator $W_h^2$ appears. The length-exchange argument explains this exceptional invariant, and the group analysis proves that the generators are complete.

At $h=0.1$, the classification shows that no independent common polynomial invariant adds information beyond $J_h$. We also established the exact identity $E_n-E_0=C_0-C_n$, connecting modified-energy conservation to physical-energy error. These results identify both the conservation supplied by the update and the correction whose variation controls the error.

The original CVT energy observation remains open. Section 5.2 identifies the driver-dependent estimate needed to explain it at $h=0.1$.