# Doubts on the solution of a differential equation

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

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 = x0
object.position = y0
object.position = 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 = x1 # X-axis
object.position = y1 # Y-axis
object.position = 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 = x0
object.position = 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 = x1 # X-axis
object.position = y1 # Y-axis
$$$$

• So you have an object affixed to a spiral path and under gravity, and you want to compute its motion by computing the force along the tangent of the spiral? Dec 13, 2022 at 8:08
• Yes, exactly. That's another way to look at it. Dec 13, 2022 at 13:54
• Please change your $t$s to $\theta$s, or vice versa. Dec 16, 2022 at 0:22
• @TedShifrin Done, I hadn't realized about that mistake in the inital $\vec{r}$. Thanks for pointing it out! Dec 17, 2022 at 12:48

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.

• Your implementation is not compatible with the cylinder coordinates in the formula. Your initial $x,y,z$ are not on the spiral of the formulas. // The stable equilibrium of the angle equation in the code is at $θ=0$ plus periods, so $x=-R\sinθ$, $z=R-R\cosθ$ makes sense with that and starting with $x=y=z=0$. But then $x(0)=0$. // With $θ(0)=0=\dotθ(0)$ the system is at rest, so nothing should happen in the solution. // Note that I changed the effective length in the formula. Dec 17, 2022 at 15:25
• You are right!! I have now adapted the code to make sure that the initial positions are on the spiral. That was a silly mistake! // As I considered for $\theta=0$ to be horizontal, I made sure that the object started at the top with a dtheta` not equal to zero. I also used the formulas for the position with respect to the angle directly instead of using the respective velocities with Euler's method. Now everything works perfectly well! Thank you for your help! (I hope you don't mind that I add you to the aknowledgments from my research paper) Dec 17, 2022 at 20:37
• Hi Lutz, I was looking through your explanation, and I'm not understanding how you can say that $(\vec{q}''(\theta)\cdot\vec{q}'(\theta))=0$ (by the way, what is $\vec{q}$? is it the $\vec{r}$ I defined above, or only a component of it?). Thanks in advance Dec 18, 2022 at 9:21
• Thanks, now I get it. $(\vec{v}\cdot\vec{q}'(\theta)) = (0, 0, g)\cdot(-R\sin, k, R\cos) = gR\cos$ (considering that $g$ is already negative) Dec 18, 2022 at 12:34
• You need the chain rule. The path is described by $x(θ)$. You want to integrate to get $X(t)=x(θ(t))$, so that $\dot X(t)=x'(θ(t)){\bf \dot θ(t)}$. This you need also to apply in the first example. Dec 18, 2022 at 14:41