I was interested in demonstrating the notion of geodesics in constrained motion and prepared the calculation of force-free motion on the unit sphere, following Hertz' take on mechanics. I wanted to set up the system in x,y,z rather than constraining motion to $r=1$ in a spherical coordinate system.

To do that we consider the movement of a point mass constrained to the unit sphere in $\mathbb{R}^3$. We let the point be a free mass located at $\boldsymbol{r}(x,y,z)=x(t)\mathbf{e}_x+y(t)\mathbf{e}_y+z(t)\mathbf{e}_z$ with the holonomic constraint $$g(x,y,z)=x^2+y^2+z^2-1.$$

We get then

\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial L}{\partial \dot{\mathbf{r}}}\right) - \frac{\partial L}{\partial \mathbf{r}} &= \lambda(t)\nabla g(x,y,z) \\ m\ddot{\mathbf{r}} &= \lambda(t) \nabla g(x,y,z) \end{align}

which delivers 4 equations

\begin{align} m\ddot{x} -2\lambda x &=0 \\ m\ddot y -2 \lambda y &=0 \\ m\ddot z -2 \lambda z &=0 \\ x^2+y^2+z^2-1&=0. \end{align}

We multiply each equation by its coordinate and add together \begin{align} m\ddot{x}x+m\ddot yy+m\ddot zz &= 2\lambda (x^2+y^2+z^2) \\ \frac m2(\ddot{x}x+\ddot yy+\ddot zz ) &= \lambda. \end{align}

Using the second time derivative of $g$ we have $$\dot x^2+\dot y^2 +\dot z^2 = -(\ddot{x}x+\ddot yy+\ddot zz ) $$ we obtain $$\lambda = -\frac m2(\dot x^2+\dot y^2 +\dot z^2) = -\frac m2 \|\dot{\boldsymbol{r}}\|^2 $$ to get \begin{align} m\ddot{x} +m(\dot x^2+\dot y^2 +\dot z^2) x &=0 \\ m\ddot y +m(\dot x^2+\dot y^2 +\dot z^2) y &=0 \\ m\ddot z + m(\dot x^2+\dot y^2 +\dot z^2) z &=0 \\ \end{align}

or $$m\ddot{\mathbf{r}} +m\|\dot{\mathbf{r}}\|^2 \mathbf{r} = 0. $$

To demonstrate I integrated this first with a leapfrog integrator and found my point mass was not staying on the sphere. After calling myself completely into question, I decided to try the Yoshida 4th-order integrator and got near-perfect results (constraints, conservation, etc.). I tried the problem in spherical coords and was also not conserving energy (although constraints satisfied).

What is so high-order about the problem that second-order integration won't work? How would one analyse the needed integration order?

I was interested in demonstrating the notion of geodesics in constrained motion and prepared the calculation of force-free motion on the unit sphere, following Hertz' take on mechanics. Since the constraints could be arbitrary and the principal would still apply, I wanted to set up the system in x,y,z rather than constraining motion to $r=1$ in a spherical coordinate system.

To do that we consider the movement of a point mass constrained to the unit sphere in $\mathbb{R}^3$. We let the point be a free mass located at $\boldsymbol{r}(x,y,z)=x(t)\mathbf{e}_x+y(t)\mathbf{e}_y+z(t)\mathbf{e}_z$ with the holonomic constraint $$g(x,y,z)=x^2+y^2+z^2-1.$$

We get then

\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial L}{\partial \dot{\mathbf{r}}}\right) - \frac{\partial L}{\partial \mathbf{r}} &= \lambda(t)\nabla g(x,y,z) \\ m\ddot{\mathbf{r}} &= \lambda(t) \nabla g(x,y,z) \end{align}

which delivers 4 equations

\begin{align} m\ddot{x} -2\lambda x &=0 \\ m\ddot y -2 \lambda y &=0 \\ m\ddot z -2 \lambda z &=0 \\ x^2+y^2+z^2-1&=0. \end{align}

We multiply each equation by its coordinate and add together \begin{align} m\ddot{x}x+m\ddot yy+m\ddot zz &= 2\lambda (x^2+y^2+z^2) \\ \frac m2(\ddot{x}x+\ddot yy+\ddot zz ) &= \lambda. \end{align}

Using the second time derivative of $g$ we have $$\dot x^2+\dot y^2 +\dot z^2 = -(\ddot{x}x+\ddot yy+\ddot zz ) $$ we obtain $$\lambda = -\frac m2(\dot x^2+\dot y^2 +\dot z^2) = -\frac m2 \|\dot{\boldsymbol{r}}\|^2 $$ to get \begin{align} m\ddot{x} +m(\dot x^2+\dot y^2 +\dot z^2) x &=0 \\ m\ddot y +m(\dot x^2+\dot y^2 +\dot z^2) y &=0 \\ m\ddot z + m(\dot x^2+\dot y^2 +\dot z^2) z &=0 \\ \end{align}

or $$m\ddot{\mathbf{r}} +m\|\dot{\mathbf{r}}\|^2 \mathbf{r} = 0. $$

To demonstrate I integrated this first with a leapfrog integrator and found my point mass was not staying on the sphere, but was converging to a wrong result. After calling myself completely into question, I decided to try the Yoshida 4th-order integrator and got near-perfect results (constraints, conservation, etc.). I tried the problem in spherical coords and was also not conserving energy (although constraints satisfied).

What is so high-order about the problem that second-order integration won't work? How would one analyse the needed integration order?

I was interested in demonstrating the notion of geodesics in constrained motion and prepared the calculation of force-free motion on the unit sphere, following Hertz' take on mechanics. I wanted to set up the system in x,y,z rather than constraining motion to $r=1$ in a spherical coordinate system.

To do that we consider the movement of a point mass constrained to the unit sphere in $\mathbb{R}^3$. We let the point be a free mass located at $\boldsymbol{r}(x,y,z)=x(t)\boldsymbol{e}_x+y(t)\boldsymbol{e}_y+z(t)\boldsymbol{e}_z$ with the holonomic constraint $$g(x,y,z)=x^2+y^2+z^2-1.$$

We get then

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial L}{\partial \dot{\boldsymbol{r}}}\right) - \frac{\partial L}{\partial \boldsymbol{r}} = \lambda(t)\nabla g(x,y,z)$$ $$ m\ddot{\boldsymbol{r}} = \lambda(t) \nabla g(x,y,z) $$

which delivers 4 equations

$$ m\ddot{x} -2\lambda x =0 $$ $$m\ddot y -2 \lambda y =0 $$ $$m\ddot z -2 \lambda z =0 $$ $$x^2+y^2+z^2-1=0. $$

We multiply each equation by its coordinate and add together $$\begin{align} m\ddot{x}x+m\ddot yy+m\ddot zz &= 2\lambda (x^2+y^2+z^2) \\ \frac m2(\ddot{x}x+\ddot yy+\ddot zz ) &= \lambda. \end{align} $$

Using the second time derivative of $g$ we have $$\dot x^2+\dot y^2 +\dot z^2 = -(\ddot{x}x+\ddot yy+\ddot zz ) $$ we obtain $$\lambda = -\frac m2(\dot x^2+\dot y^2 +\dot z^2) = -\frac m2 \|\dot{\boldsymbol{r}}\|^2 $$ to get $$\begin{align} m\ddot{x} +m(\dot x^2+\dot y^2 +\dot z^2) x &=0 \\ m\ddot y +m(\dot x^2+\dot y^2 +\dot z^2) y &=0 \\ m\ddot z + m(\dot x^2+\dot y^2 +\dot z^2) z &=0 \\ \end{align} $$ or $$m\ddot{\boldsymbol{r}} +m\|\dot{\boldsymbol{r}}\|^2 \boldsymbol{r} = 0. $$

To demonstrate I integrated this first with a leapfrog integrator and found my point mass was not staying on the sphere. After calling myself completely into question, I decided to try the Yoshida 4th-order integrator and got near-perfect results (constraints, conservation, etc.). I tried the problem in spherical coords and was also not conserving energy (although constraints satisfied).

What is so high-order about the problem that second-order integration won't work? How would one analyse the needed integration order?

I was interested in demonstrating the notion of geodesics in constrained motion and prepared the calculation of force-free motion on the unit sphere, following Hertz' take on mechanics. I wanted to set up the system in x,y,z rather than constraining motion to $r=1$ in a spherical coordinate system.

To do that we consider the movement of a point mass constrained to the unit sphere in $\mathbb{R}^3$. We let the point be a free mass located at $\boldsymbol{r}(x,y,z)=x(t)\mathbf{e}_x+y(t)\mathbf{e}_y+z(t)\mathbf{e}_z$ with the holonomic constraint $$g(x,y,z)=x^2+y^2+z^2-1.$$

We get then

\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial L}{\partial \dot{\mathbf{r}}}\right) - \frac{\partial L}{\partial \mathbf{r}} &= \lambda(t)\nabla g(x,y,z) \\ m\ddot{\mathbf{r}} &= \lambda(t) \nabla g(x,y,z) \end{align}

which delivers 4 equations

\begin{align} m\ddot{x} -2\lambda x &=0 \\ m\ddot y -2 \lambda y &=0 \\ m\ddot z -2 \lambda z &=0 \\ x^2+y^2+z^2-1&=0. \end{align}

We multiply each equation by its coordinate and add together \begin{align} m\ddot{x}x+m\ddot yy+m\ddot zz &= 2\lambda (x^2+y^2+z^2) \\ \frac m2(\ddot{x}x+\ddot yy+\ddot zz ) &= \lambda. \end{align}

Using the second time derivative of $g$ we have $$\dot x^2+\dot y^2 +\dot z^2 = -(\ddot{x}x+\ddot yy+\ddot zz ) $$ we obtain $$\lambda = -\frac m2(\dot x^2+\dot y^2 +\dot z^2) = -\frac m2 \|\dot{\boldsymbol{r}}\|^2 $$ to get \begin{align} m\ddot{x} +m(\dot x^2+\dot y^2 +\dot z^2) x &=0 \\ m\ddot y +m(\dot x^2+\dot y^2 +\dot z^2) y &=0 \\ m\ddot z + m(\dot x^2+\dot y^2 +\dot z^2) z &=0 \\ \end{align}

or $$m\ddot{\mathbf{r}} +m\|\dot{\mathbf{r}}\|^2 \mathbf{r} = 0. $$

To demonstrate I integrated this first with a leapfrog integrator and found my point mass was not staying on the sphere. After calling myself completely into question, I decided to try the Yoshida 4th-order integrator and got near-perfect results (constraints, conservation, etc.). I tried the problem in spherical coords and was also not conserving energy (although constraints satisfied).

What is so high-order about the problem that second-order integration won't work? How would one analyse the needed integration order?