computationalphysics

chaotic tumbling of hyperion

abstract

as wahat we already know pyperion,one smaller moons of Saturn,is tumbling in a chaotic way due to its irregular shape.I will study the simplified hyperion model in a nimerical way,and estimate the lyapunov exponenent for divergence of two nearby trajectories of the tumbling motion of of Hyperion.

background

Hyperion, also known as Saturn VII, is a moon of Saturn discovered by William Cranch Bond, George Phillips Bond and William Lassell in 1848. It is distinguished by its irregular shape, its chaotic rotation, and its unexplained sponge-like appearance. It was the first non-round moon to be discovered.

- equations for Hyperion's trajectory
the inverse-square law of gravitation:

$$\vec{F}=G\frac{Mm}{r^3}\hat{r}$$
we write it as two first-order differental equation so that we can use euler-cromer method: $$\frac{dv_x}{dt}=-\frac{GM_{sat} x}{r_c^3}\\ \frac{dx}{dt}=v_x\\ \frac{dv_y}{dt}=-\frac{GM_{sat} y}{r^3}\\ \frac{dy}{dt}=v_y$$
for circular motion we know that the force must be equal to $M_P v^2 /r$,which lead to:
$$\frac{M_H v^2}{r}=F_G=\frac{GM_{sat} M_{H}}{r^2}$$
where we chose proper units(unit of lenth to be radius of Hyperion's orbit,which we might called $1\space HU=$"Hyperion unit",and that of time to be the orbit period of Hyperion around Saturn(1"Hyperion-year")):
$$GM_{sat}=v_h^2r_h=4\pi r^2\space HU^3/Hr^2$$
Given the equations above and Euler-cromer method, we can calculate the trajectory of Hyperion.
- equations for qualitively estimate the chaotic tumbling of Hyperion
since our objective is simply to show that teh motion of such an irregularly shaped moon can be chaotic,we consider the model as two particles ,$m_1$ and $m_2$,connected bby a massless rod in orbit around a massive object located at the origion..we let $\theta$ be the angle that the rod makes with the $x$ axis and define the associated angular velocity $\omega=\frac{d\theta}{dt}$,and $r_c$ as the distence of the motion of the center of the mass ,$x_c$ and $y_c$ as the position of teh center of mass.
after some algebra,we get: $$\frac{d\omega}{dt}\approx-\frac{3GM_{sat}}{r_c^5}(x_csin\theta-y_ccos\theta)(x_ccos\theta+y_csin\theta)$$

main body

• the python programme

import numpy as np
import pylab as pl
import math
class hyperion:
    def __init__(self,x0=1,y0=0,vx=0,vy=2*math.pi,dt=0.0001,total_time=10,initial_theta=0,initial_omega=0):
        self.x=[x0]
        self.y=[y0]
        self.r=[math.sqrt(x0**2+y0**2)]
        self.vx=[vx]
        self.vy=[vy]
        self.dt=dt
        self.t=[0]
        self.tt=total_time
        self.th=[initial_theta]
        self.om=[initial_omega]
        self.a=0
        self.b=0
        self.ecc=0
        self.dtheta=[]
    def run(self):
        while self.t[-1]<self.tt:
            self.vx.append(self.vx[-1]-4*math.pi**2*self.x[-1]/self.r[-1]**3*self.dt)
            self.vy.append(self.vy[-1]-4*math.pi**2*self.y[-1]/self.r[-1]**3*self.dt)
            self.x.append(self.x[-1]+self.vx[-1]*self.dt)
            self.y.append(self.y[-1]+self.vy[-1]*self.dt)
            self.r.append(math.sqrt(self.x[-1]**2+self.y[-1]**2))
            temp=-12*math.pi**2*(self.x[-1]*math.sin(self.th[-1])-self.y[-1]*math.cos(self.th[-1]))\
            *(self.x[-1]*math.cos(self.th[-1])+self.y[-1]*math.sin(self.th[-1]))
            self.om.append(self.om[-1]+temp*self.dt)
            self.th.append(self.th[-1]+self.om[-1]*self.dt)
            if self.th[-1]>math.pi:
                self.th[-1]-=2*math.pi
            elif self.th[-1]<-math.pi:
                self.th[-1]+=2*math.pi
            self.t.append(self.t[-1]+self.dt)
        self.a=(max(self.x)-min(self.x))/2
        self.b=(max(self.y)-min(self.y))/2
        self.ecc=math.sqrt(math.fabs(self.a**2-self.b**2))/self.a
    def showth(self):
        pl.title("Hyperion $\\theta$ versus time")
        pl.xlabel("time(Hyperion-year)")
        pl.ylabel("$\\theta (radius)$")
        pl.xlim(0,10)
        pl.plot(self.t, self.th,label="eccentricity=%.2f"%self.ecc) 
        pl.legend()
        pl.grid(True)
        pl.show()
    def showom(self):
        pl.title("Hyperion $\\omega$ versus time")
        pl.xlabel("time(Hyperion-year)")
        pl.ylabel("$\\omega (radius)$")
        pl.xlim(0,10)
        pl.plot(self.t, self.om,label="eccentricity=%.2f"%self.ecc)  
        pl.legend()
        pl.grid(True)
        pl.show()
#plot theta or omega versus time
a=hyperion()
a.run()
a.showth()
a=hyperion()
a.run()
a.showom()
#plot delta theta versus time and calculate the lyapuniv exponent
class qwer(hyperion):
    def haha(self):
        _d=0.01
        _vy=2*math.pi-1.1
        a=hyperion(vy=_vy)
        a.run()
        print(a.ecc)
        tempth1=a.th
        a=hyperion(vy=_vy,initial_theta=_d)
        a.run()
        tempth2=a.th
        for i in range(len(tempth2)):
            a.dtheta.append(math.log(math.fabs(tempth2[i]-tempth1[i])))
        pl.plot(a.t,a.dtheta,label="eccentricity=%.4f"%a.ecc)
        dd=[]
        tt=[]
        for j in range(len(a.dtheta)):
            dd.append(a.dtheta[j])
            tt.append(a.t[j])
        z=np.polyfit(tt,dd,1)
        p=np.poly1d(z)
        print(p)
        linspx=np.linspace(0,8)
        linspy=z[0]*linspx+z[1]
        pl.plot(linspx,linspy,label="lyapunov exp=%.4f"%z[0])
        pl.xlabel("time(Hyperion-year)")
        pl.ylabel("$\\Delta\\theta(radius)$")
        pl.title("Hyperion $\\Delta\\theta$ versus time for initial $\\Delta\\theta$=%.3f"%_d)
        pl.legend()
b=qwer()
b.haha()

• 1.Using this python programme,we can conclude that,as the eccentricity is not zero,the system become more and more chaotic.
  
• How the lyapunov exponnet varies with the eccentricity
  In the plots above,we can see that when the eccentricity=0.00(actually its 0.0004,you can see this if you run my programme by yourself),we barely see any chaotic phenomenon in Hyperion's tumbling.but when we come to calculate the lyapunov exponent,we can obviously see the lyapunov exponent is not zero but 0.02.
  

eccentricity	lyapunov exponent		eccentricity	lyapunov exponent
0.0004	0.0234		0.1528	0.4547
0.0316	0.0229		0.1819	0.5106
0.0626	0.0133		0.2104	0.3092
0.0932	0.0129		0.2384	0.2769
0.1223	0.0699		0.2930	0.1793

I plot the table in the next picture:

it doesn't looks like anything to me.

• 3.another important factor that can not be ignored.
  This factor is the step-time we used .When we choose the units for lenth is $1HU$,when the initial velocity is $2\pi$,according to theoretical analysis,the trajectory is suppose to be a circular orbit and Hyperion tumbling in a regular way,which means the eccentricity is zero.But the programme turns out to be 0.0004 and Hyperion tumbling in a chaotic way.so we can conclude that the step-time have significiant influence on the simulation results,which tell us a smaller step-time is neccessary.

conclusion

• 1.A moon which have a irregular shape most probably tumbling chaoticly.
• 2.The step-time have significiant influence on the simulation results,which tell us a smaller step-time is neccessary.

acknowledgement

Nicholas J.Giodano's computational physics.