Starfish Orbit

This was fun! Presenting my numerical solution.

Consider the system at a general instant of time. The position of the particle is:

$x = r\cos{\theta}$ $y = r\sin{\theta}$

$\dot{x} = \dot{r}\cos{\theta} - r\dot{\theta}\sin{\theta}$ $\dot{y} = \dot{r}\sin{\theta} + r\dot{\theta}\cos{\theta}$

The kinetic energy of the system is:

$T = \frac{1}{2}\left(\dot{x}^2 + \dot{y}^2\right) = \frac{1}{2}\left(\dot{r}^2 + r^2 \dot{\theta}^2\right)$ $V = -r^{-1} + r^{-3}$

Applying Lagrange's equation for the $\theta$ coordinate gives:

$\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{\theta}}\right) - \frac{\partial T}{\partial \theta} + \frac{\partial V}{\partial \theta}=0$ $\implies \frac{d}{dt}\left(r^2 \dot{\theta}\right) = 0$ $\implies r^2\dot{\theta} = K \ \dots \ (1)$

Where $K$ is some constant. This constant can be found by applying the initial conditions which gives:

$\dot{\theta} = \frac{2}{5r^2} \ \dots (2)$

Applying Lagrange's equation for the $r$ coordinate gives: $\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{r}}\right) - \frac{\partial T}{\partial r} + \frac{\partial V}{\partial r}=0$ $\ddot{r} = r\dot{\theta}^2 -\frac{1}{r^2} + \frac{3}{r^4}$

Replacing (2) in the equation above leads to a differential equation:

$\ddot{r} = \frac{4}{25 r^3} -\frac{1}{r^2} + \frac{3}{r^4} \ \dots (3)$

So finally, the system is governed by the differential equations:

$\boxed{\dot{\theta} = \frac{2}{5r^2} }$ $\boxed{\ddot{r} = \frac{4}{25 r^3} -\frac{1}{r^2} + \frac{3}{r^4} }$

Using the given initial conditions, the above equations are solved numerically. First, the system is simulated until $t=60$ . The resulting trajectory is:

The trajectory above indeed looks like a ruptured starfish.

Simulating till $t=120$ yields:

One sees that the trajectory looks like two ruptured starfish one atop another. It gets really interesting when the system is simulated until $t=1000$

So, what one observes here is that the values of $r$ which the system solution obtains, are confined to a lower and upper bound. $r$ varies between $1.501$ and $2.3356$ . If one were to look at the plot of time vs. the x or y coordinate, it can be seen that the motion shows a semblance of periodic repeatability, but not exact periodicity. I think that the technical term for this is quasi-periodicity.

A very well known quasi-periodic dynamic system to us is the earth's climate and how it exhibits yearly repeatability which we understand as seasons, but the periodicity is not exact. There are always slight variations in the time of onset and intensity of various seasons.

If one were to simulate the given system for a significantly longer duration (say $t=20000$ ), it is seen that every point within the bounded area would be reached by the system. The trajectory would look like a perfect blue doughnut as shown below:

As for why this system exhibits such behaviour, I have an explanation, which I am not sure is correct. Consider the equation of motion:

$\ddot{r} = \frac{4}{25 r^3} -\frac{1}{r^2} + \frac{3}{r^4}$

Here $\ddot{r} =0$ when $r \approx 1.813$ . So about this equilibrium, which happens to be a stable equilibrium state for the $r(t)$ solution, the system oscillates in an almost periodic manner, or a quasi-periodic manner.

I use Mathcad.

Next fine problems -

$1.$ find initial data for choreograph trajectories (examples here) $2.$ find orbit for initial data

$t=0 \\ r=\displaystyle\frac{2+\sqrt{1879}}{25} \\ \dot{\theta}=\frac{250}{(2+\sqrt{1879})^{2}}$

$3.$ find the analitic equation for $r_{min}$ and $r_{max}$ (for this problem r varies between $r_{min} =1.501$ and $r_{max}=2.3356$ ).

As some background, I initially had two attractive potentials, but the particle kept spiraling in toward the origin and causing numerical problems. So I then made it so that the potential is attractive for "large" radius and repulsive for "small" radius. That dual nature of the potential is what creates the starfish pattern.

@Karan Chatrath has already provided an excellent solution based on Lagrangian mechanics, which is how I initially did it. I also did it the Newtonian way, just for fun. Initialize the $(\dot{x}, \dot{y})$ values based on the $(\dot{r}, \dot{\theta})$ at the beginning:

At time $t = 0$ :

$\dot{x} = \dot{r} \cos \theta - r \sin \theta \, \dot{\theta} \\ \dot{y} = \dot{r} \sin \theta + r \cos \theta \, \dot{\theta}$

Then we can derive the net force on the particle as follows:

$\vec{v} = (x,y) \\ \vec{u} = \frac{\vec{v}}{|\vec{v}|} = \frac{\vec{v}}{r} \\ F = -\frac{dV}{dr} = -r^{-2} + 3 r^{-4} \\ \vec{F} = F \vec{u}$

And then we can numerically integrate as usual, and the result is identical to the one derived using Lagrange, as expected. The Newton solution code is attached.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

import math

# constants

m = 1.0
dt = 10.0**(-5.0)

##################################################

# Initial conditions in (r,theta)

t = 0.0
count = 0

r = 2.0
theta = 0.0

rd = 0.2
thetad = 0.1

##################################################

# Initial conditions in (x,y)

x = r*math.cos(theta)
y = r*math.sin(theta)

xd = rd*math.cos(theta) - r*math.sin(theta)*thetad
yd = rd*math.sin(theta) + r*math.cos(theta)*thetad

xdd = 0.0
ydd = 0.0

##################################################

while t <= 60.0:

    x = x + xd*dt    # numerical integration in (x,y)
    y = y + yd*dt

    xd = xd + xdd*dt
    yd = yd + ydd*dt

    r = math.hypot(x,y)

    ux = x/r
    uy = y/r

    F = -(r**-2.0) + 3.0*(r**-4.0)

    Fx = F*ux
    Fy = F*uy

    xdd = Fx/m
    ydd = Fy/m

    t = t + dt
    count = count + 1

    if count % 1000 == 0:

        print t,x,y,r

##################################################