Calculate u in terms of time such that a particle maintains a constant speed following a parametric equation

Question

I have a parametric equation given by:

$x=\cos(6u)$
$y=\sin(4u)$

And I understand that the speed of a particle at any given t is:

$\sqrt{\left(\dfrac{dx}{du}\right)^2 + \left(\dfrac {dy}{du} \right)^2} $

Which for my parametric equation would be:

$\sqrt{36 \sin^2(6u) + 16 \cos^2(4u)}$

I am trying to make a particle move at a constant speed along the parametric equation.

I have looked and seem to be unable to find information on how to do this, the only thing I did find is this paper on the subject.

Unfortunately given my rather basic grasp of calculus I wasn't able to make much sense of it.

To clarify, I am looking for $u$ in terms of time such that the particle maintains constant speed.

Thanks.

Answer 1

Generalities: Suppose a plane parametric curve is defined by $x = f(t)$ and $y = g(t)$ for $a \leq t \leq b$. (For technical reasons, assume $f$ and $g$ have continuous first derivatives, and that the curve's velocity is never $(0, 0)$, i.e., $f'$ and $g'$ never vanish at the same point.)

The speed of the curve is $$ s'(t) = \sqrt{f'(t)^{2} + g'(t)^{2}}, $$ and the arc length function $s$ on $[a, b]$ is defined by $$ s(t) = \int_{a}^{t} \sqrt{f'(u)^{2} + g'(u)^{2}}\, du. $$ For convenience, let $\ell = s(b)$ denote the length of the curve.

Since $s' > 0$ (the velocity is non-vanishing, so the speed is positive), $s:[a, b] \to [0, \ell]$ is strictly increasing, and therefore has a differentiable inverse function $\tau:[0, \ell] \to [a, b]$. By the chain rule and the formula for the derivative of an inverse function, the parametrization $$ x = f \circ \tau,\quad y = g \circ \tau $$ has speed $s'(t)\, \tau'(s) = 1$.

---

Specifics: For the curve in this question, $x = \cos(6t)$ and $y = \sin(4t)$ are $\pi$-periodic, and the arc length function is given by $$ s(t) = \int_{0}^{t} \sqrt{36 \sin^{2}(6u) + 16 \cos^{2}(4u)}\, du,\quad 0 \leq t \leq \pi. $$ As above, put $$ \ell = s(\pi) = \int_{0}^{\pi} \sqrt{36 \sin^{2}(6u) + 16 \cos^{2}(4u)}\, du. $$

Let's suppose the animation is to have $N$ "position steps", i.e., the particle can be at any of $N$ distinct locations. (In the diagram, $N = 120$.)

One (inflexible, and therefore possibly inelegant) approach would be to pre-compute (and store in the program) $N$ values $t_{i}$, $1 \leq i \leq N$, for which $$ s(t_{i}) \approx \frac{i}{N}\, \ell, $$ namely, to calculate the $t$ values for which the particle is $(i/N)$ of the way along the curve with respect to arc length. Then, when the program needs to display the particle at "the $i$th position", it plots $\bigl(f(t_{i}), g(t_{i})\bigr)$. (Symmetry of the speed for this curve allows the position to be reconstructed from $N/4$ points, in case that's helpful.)

---

Note: The original parametrization, the "shadow" of a point moving at constant speed along an imaginary surface (a flat torus in four-dimensional space if it matters), may look more "natural" than a constant-speed parametrization. That is, don't expend a lot of computational effort to make the motion look "forced" unless you're positive that's what you want or need.

Answer 2

An easy numerical solution is obtained by solving the differential equation for $u$ as a function of $t$:

$$V=\frac{ds}{dt}=\frac{ds}{du}\frac{du}{dt},$$ or $$\frac{du}{dt}=\frac V{\sqrt{36 \sin^2(6u) + 16 \cos^2(4u)}}.$$

Use the Runge-Kutta method.

Answer 3

As Yves Daoust explained, you can use one of the Runge-Kutta methods to numerically solve the differential equation for your parameter. However, if the step size is small (the speed is low) you can get by with simple Euler integration. In other words,

$\frac{\Delta u}{\Delta t} \approx \frac{du}{dt}$

so just let $\Delta u = \Delta t \cdot \frac{du}{dt}$

To optimize this calculation we can express the arc length differential in terms of $x$ and $y$. The equations below use $a$ and $b$ for the $x$ & $y$ frequencies, rather than your 6 & 4, in order to make the result apply to any Lissajous figure. They also use $\theta$ for the parameter; I figured that was appropriate because we're dealing with trigonometric functions, and I guess it looks a little more professional. :)

$$\begin{align} \text{Let } x & = \cos (a \theta), \, y = \sin (b \theta) \\ \frac{dx}{d\theta} & = - a \sin (a \theta), \, \frac{dy}{d\theta} = b \cos (b \theta) \\ \\ \frac{ds}{d\theta}^2 & = \frac{dx}{d\theta}^2 + \frac{dy}{d\theta}^2 \\ & = a^2 \sin^2 (a \theta) + b^2 \cos^2 (b \theta) \\ & = a^2 (1 - \cos^2 (a \theta)) + b^2 (1 - \sin^2 (b \theta)) \\ & = a^2 (1 - x^2) + b^2 (1 - y^2) \\ \frac{ds}{d\theta} & = \sqrt{a^2 (1 - x^2) + b^2 (1 - y^2)}\\ \\ v & = \frac{ds}{dt} = \frac{ds}{d\theta} \cdot \frac{d\theta}{dt}\\ \\ \text{Therefore}\\ \frac{d\theta}{dt} & = \frac{v}{ds/d\theta}\\ \Delta\theta & \approx \frac{\Delta t \cdot v}{ds/d\theta}\\ & = \frac{\Delta t \cdot v}{\sqrt{a^2 (1 - x^2) + b^2 (1 - y^2)}} \end{align}$$

Here's a small Python snippet to illustrate the algorithm:

def liss(a, b, speed, width, height):
    pix = bitmap(width, height)
    ox, oy = width // 2, height // 2
    rad = min(ox, oy) - 3

    #Precompute constants used in speed calculation
    r2 = rad * rad
    a2 = a * a
    b2 = b * b
    c2 = (a2 + b2) * r2

    maxtheta = 2. * pi
    theta = 0.0
    while theta < maxtheta:
        x = rad * cos(a * theta)
        y = rad * sin(b * theta)
        pix.put(int(.5 + ox + x), int(.5 + oy + y))
        dtheta = speed / (c2 - a2 * x*x - b2 * y*y) ** 0.5
        theta += dtheta

Here's the output of the above algorithm, with a=3, b=2, speed=10:

As you can see, it's fairly smooth. But we can make it better by following Yves' advice and using Runge-Kutta integration. Here's Python code that uses RK4:

def liss(a, b, speed, width, height):
    pix = bitmap(width, height)
    ox, oy = width // 2, height // 2
    rad = min(ox, oy) - 3

    def speedfunc(theta):
        s = (a * sin(a * theta))**2 + (b * cos(b * theta))**2
        return speed / (rad * s ** 0.5)

    maxtheta = 2. * pi
    theta = 0.0
    while theta < maxtheta:
        x = rad * cos(a * theta)
        y = rad * sin(b * theta)
        pix.put(int(.5 + ox + x), int(.5 + oy + y))

        #Calculate Runge-Kutta 4 increments
        k1 = speedfunc(theta)
        k2 = speedfunc(theta + 0.5 * k1)
        k3 = speedfunc(theta + 0.5 * k2)
        k4 = speedfunc(theta + k3)

        dtheta = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
        theta += dtheta

And its output, with the same parameters as above:

It's hard to see much difference here, but the RK4 algorithm is clearly superior for high speeds, especially in the regions of the curve where the acceleration is high.

In case you're interested, here's a Pastebin link to the full Python 2 program used to produce the RK4 version as a plain PBM file; you can easily adapt it to use the Euler algorithm. This program was tested on Python 2.6 but it probably runs on earlier versions; it will need minor modifications to run on Python 3.

LissRK4.py

Finally, I'd like to mention that although these constant speed algorithms make nice-looking plots, I have to agree with user86418 that the simple constant angle increment algorithm looks more natural in an animation; after all, this is simple harmonic motion. :)