Problem 4.19：Study the behavior of our model for Hyperion for different initial conditions.Estimate the Lyapunov exponent .Examine how this exponent varies as a function of the eccentricity of the orbit.
background:
The Main Body:
The model has been simplified as figure 4.16 in textbook.There are two forces acting on each of the masses, the force of the gravity from Saturn and the force from the rod. Since we are interested in the motion about the center of mass, the force from the rod does not contribute. The gravitational force on m1 can be written as
where Msat is the mass of Saturn, r1 is the distance from Saturn to m1, i and j are unit vectors in the x and y directions. The coordinates of the center of mass are (xc,yc), so that (x1-xc)i+(y1-yc)j is the vector from the center of mass to m1. The torque on m1 is then
with a similar expression for r2. The total torque on the moon is just r1+r2 and this is related to the time derivative of omega by:
Hence,we have
rc is the distance from the center of mass to Saturn.
Code:
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()
as the eccentricity is not zero,the system become more and more chaotic:
How the lyapunov exponnet varies with the eccentricity:
Conclusion:
This extreme sensitivity to initial conditions is one of the hallmarks of chaotic behavior.
