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
```
Best Answer
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.