Doubts on the solution of a differential equation

Question

I've encountered myself with a differential equation which I'm not sure that I've solved correctly.

The situation is the following: I have a set of parametric equations, which I'll call $\vec{r}$

\begin{array}{ll} x(\theta) &= R \cdot \cos(\theta) \cr y(\theta) &= k \cdot \theta \cr z(\theta) &= R \cdot \sin(\theta) \end{array}

$x$, $y$ and $z$ represent position in all three axes.

As I'm going to need them later, I'll calculate the first and second derivative of $\vec{r}$

\begin{array}{ll} x\prime(\theta) &= - R \cdot \sin(\theta) \cr y\prime(\theta) &= k \cr z\prime(\theta) &= R \cdot \cos(\theta) \end{array}

\begin{array}{ll} x\prime\prime(\theta) &= - R \cdot \cos(\theta) \cr y\prime\prime(\theta) &= 0 \cr z\prime\prime(\theta) &= - R \cdot \sin(\theta) \end{array}

The resulting directional vector would be written as follows: $\vec{u} = (- R \cdot \sin(\theta), k, R \cdot \cos(\theta))$. Using the Pythagorean theorem, we can calculate the modulus of this vector: $\|\vec{u}\| = \sqrt{R^2 \cdot \sin^2(\theta) + k^2 + R^2 \cdot \cos^2(\theta)}$. Combining these, we obtain the unitary vector:

$$\hat{u} = \frac{\vec{u}}{\|\vec{u}\|} = \frac{(- R \cdot \sin(\theta), k, R \cdot \cos(\theta))}{\sqrt{R^2 + k^2}}$$

Therefore, using the vector for the gravitational acceleration $\vec{v} = (0, 0, g)$ and knowing that $\|\vec{a}_T\| = \|\vec{v}\| \cdot \cos(\beta) = \|(\vec{v} \cdot \hat{u})\|$, we can obtain the modulus of the total acceleration:

$$\|\vec{a}_T\| = \frac{\sqrt{(0^2 + 0^2 + g^2 \cdot R^2 \cdot \cos^2(\theta))}}{\sqrt{R^2 + k^2}} = \frac{g \cdot R \cdot \cos(\theta)}{\sqrt{R^2 + k^2}}$$

The total acceleration will then result in multiplying this modulus we have just obtained by the unit vector which points towards the direction of the acceleration:

$$\vec{a} = \|\vec{a_T}\| \cdot \hat{u} = \frac{g \cdot R \cdot \cos(\theta)}{\sqrt{R^2 + k^2}} \cdot \frac{(- R \cdot \sin(\theta), k, R \cdot \cos(\theta))}{\sqrt{R^2 + k^2}}$$

$$\vec{a} = \frac{(- g \cdot R^2 \cdot \cos(\theta) \cdot \sin(\theta), g \cdot k \cdot R \cdot \cos(\theta), g \cdot R^2 \cdot \cos^2(\theta))}{R^2 + k^2}$$

Writing it in parameterized form (I used the chain rule on the left side):

\begin{array}{rl} x: \ddot{\theta}(t) \cdot x\prime(\theta) + \dot{\theta}(t) \cdot x\prime\prime(\theta) &= \displaystyle \frac{- g \cdot \cos(\theta) \cdot \sin(\theta)}{1 + k^2} \cr y: \ddot{\theta}(t) \cdot y\prime(\theta) + \dot{\theta}(t) \cdot y\prime\prime(\theta) &= \displaystyle \frac{g \cdot k \cdot \cos(\theta) \cdot \sin(\theta)}{1 + k^2} \cr z: \ddot{\theta}(t) \cdot z\prime(\theta) + \dot{\theta}(t) \cdot z\prime\prime(\theta) &= \displaystyle \frac{g \cdot \cos^2(\theta)}{1 + k^2} \end{array}

Isolating $\ddot{\theta}(t)$ and subtituting with the formulas above, I get the following differential equations which I know how to solve numerically by using Euler's method:

\begin{array}{rl} x: \ddot{\theta}(t) &= \displaystyle \frac{- g \cdot \cos(\theta)}{R \cdot (1 + k^2)} + \dot{\theta}^2(t) \cdot \tan(\theta) \cr y: \ddot{\theta}(t) &= \displaystyle \frac{g \cdot k \cdot \cos(\theta)}{1 + k^2} \cr z: \ddot{\theta}(t) &= \displaystyle \frac{g \cdot \cos(\theta)}{R \cdot (1 + k^2)} - \dot{\theta}^2(t) \cdot \tan^{-1}(\theta) \end{array}

My question is: shouldn't the three last formulas be the same or equivalent to each other?

Edit: for anyone wondering, the final objective of this is to make an animation using Python. Here's the final code (edited out to only include the math):

import math

# VARIABLES
g = - 9.8
R = 10
k = 2
ddtheta = 0

theta0 = math.radians(90)
dtheta0 = 1
x0 = R * math.cos(theta0)
y0 = k * theta0
z0 = R * math.sin(theta0)

x1 = 0
y1 = 0
z1 = 0
theta1 = 0
dtheta1 = 0

object.position[0] = x0
object.position[1] = y0
object.position[2] = z0

# ANIMATION
for i in range(1, 200):
    t = 0.00003 # Small step size
    costheta0 = math.cos(theta0)
    sintheta0 = math.sin(theta0)
    
    # EULER'S METHOD FOR THETA
    ddtheta = (g * R * costheta0)/(R**2 + k**2)
    dtheta1 = ddtheta * t + dtheta0
    theta1 = dtheta0 * t + theta0
    
    # THE THREE AXES
    x1 = R * costheta0
    y1 = k * theta0
    z1 = R * sintheta0
        
    # RESET VALUES FOR NEXT ITERATION
    x0 = x1
    y0 = y1
    z0 = z1
    theta0 = theta1
    dtheta0 = dtheta1
    
    object.position[0] = x1 # X-axis
    object.position[1] = y1 # Y-axis
    object.position[2] = z1 # Z-axis

EDIT: The issue is solved... See Lutz's answer below...

But now, I tried applying the same method to another curve: the clothoid, but it hasn't worked. The clothoid has the following formulas for $\vec{q}$:

$$\begin{array}{rl} x(\theta) &= \displaystyle \int^\theta_0 \cos(t^2)\ dt \cr y(\theta) &= \displaystyle \int^\theta_0 \sin(t^2)\ dt \end{array}$$

The derivatives are:

$$\begin{array}{rl} x'(\theta) &= \cos(\theta^2) \cr y'(\theta) &= \sin(\theta^2) \end{array}$$

$$\begin{array}{rl} x''(\theta) &= -2 \cdot \theta \cdot \sin(\theta^2) \cr y''(\theta) &= 2 \cdot \theta \cos(\theta^2) \end{array}$$

Knowing that the chain rule is $\vec{a} = \vec{q}'' \cdot \dot{\theta}^2 + \vec{q}' \cdot \ddot{\theta}$ and that $\vec{q}'$ and $\vec{q}''$ are orthogonal ($\cos\alpha = 0$):

$$\require{cancel}\begin{array}{c} \vec{q}''(\theta) \cdot \vec{q}'(\theta) = \|\vec{q}''\| \cdot \|\vec{q}'\| \cdot \cancelto{0}{\cos 90} = 0 \cr \cr \vec{a} \cdot \vec{q}' = \underset{0}{\underbrace{\vec{q}''(\theta) \cdot \vec{q}'(\theta) \cdot \dot{\theta}^2}} + \vec{q}'(\theta)^2 \cdot \ddot{\theta}(t) \cr \cr \dfrac{\vec{F}}{m} \cdot \vec{q}'(\theta) = (0, g) \cdot (\cos(\theta^2), \sin(\theta^2)) = g \cdot \sin(\theta^2) \end{array}$$

By combining the expressions above (and knowing that the denominator $\sin^2(\theta^2) + \cos^2(\theta^2) = 1$):

$$\require{cancel}\begin{array}{c} \cancelto{1}{\|\vec{q}\|(\theta)^2} \cdot \ddot{\theta}(t) = g \cdot \sin(\theta^2) \end{array}$$

Why can't I apply the same method as explained by @Lutz Lehmann??

Here's the code for this equation:

include math
# VARIABLES
fps = 30
g = - 9.8
ddtheta = 0

theta0 = 0
dtheta0 = 5
x0 = 0
y0 = 0
vx0 = math.cos(theta0**2)
vy0 = math.sin(theta0**2)

x1 = 0
y1 = 0
vx1 = 0
vy1 = 0
theta1 = 0
dtheta1 = 0

object.position[0] = x0
object.position[2] = y0

# ANIMATION
for i in range(1, 2000):
    t = 0.00003 # Small step size
    ddtheta = g * math.sin(theta0**2)
    dtheta1 = ddtheta * t + dtheta0
    theta1 = dtheta0 * t + theta0
        
    vx1 = math.cos(theta0**2)
    x1 = vx0 * t + x0
    vy1 = math.sin(theta0**2)
    y1 = vy0 * t + y0
    
    x0 = x1
    y0 = y1
    vx0 = vx1
    vy0 = vy1
    theta0 = theta1
    dtheta0 = dtheta1
    
    object.position[0] = x1 # X-axis
    object.position[2] = y1 # Y-axis
```

Answer 1

You already laid out that the acceleration along the path has to be tangent to the path, $m\vec a=(\vec F\vcenter{\huge⋅}\hat u)\,\hat u$, as all force components perpendicular to the path are counter-acted by the path itself, be it as physical construction or force field.

As $\vec a$ is also the second time derivative of $\vec q(θ(t))$, one gets $\vec a=\vec q''(θ)\dot θ^2+\vec q'(θ)\ddot θ$, where here $(\vec q''(θ)\vcenter{\huge⋅}\vec q'(θ))=0$. For the projection we thus get $$ (\vec a \vcenter{\huge⋅}\vec q'(θ))=\|\vec q'(θ)\|^2\ddot θ \\ =(\vec F/m \vcenter{\huge⋅}\vec q'(θ))=-g·R\cos(θ) $$ Using $\|\vec q'(θ)\|^2=R^2+k^2$, this results in $$ \ddotθ= -\frac{g \cdot R \cdot \cos(\theta)}{R^2 + k^2} $$ along the tangent direction.

So what this reduces to is a pendulum equation with the increased length $\frac{R^2+k^2}R$, the spiral being oriented horizontally, so the system resembles a loop-de-loop. Note that $θ=0$ is the horizontal $x$ direction.

To get a function table in $(t,x,y,z)$ for the motion, you need to integrate this second order equation and then insert the function table for $θ(t)$ into the formula for the cylinder coordinates.

One could, as you did, add the derivatives of the cylinder coordinates to the differential system. But using the low-order Euler method would introduce additional errors, moving away from the spiral, over those from the angle integration.

As it is a Hamiltonian system, you can easily change the Euler method to semi-implicit or symplectic Euler, or with another small twist, into the leapfrog Verlet method.