Code

import math
import pylab as pl
class Hyperion:
    def __init__(self,theta1,theta2): 
        self.theta1 = [theta1]
        self.theta2 = [theta2]
        self.dt = 0.0001
        self.t = [0]
        self.x = [1]
        self.y = [0]
        self.v_x = [0]
        self.v_y = [2.0*math.pi]
        #self.v_y = [5]
        self.r = [1]
        self.beta = 2
        self.avelo1 = [0]
        self.avelo2 = [0]
        self.dthe = [math.log(abs(0.01))]
    def calculate(self):
        while (self.t[-1]<=10):
            self.r[-1] = math.sqrt((self.x[-1]**2) + (self.y[-1]**2))
            self.v_x.append(self.v_x[-1] - \
            4*(math.pi)**2*self.x[-1]*self.dt/((self.r[-1])**(self.beta + 1)))
            self.v_y.append(self.v_y[-1] - \
            4*(math.pi)**2*self.y[-1]*self.dt/((self.r[-1])**(self.beta + 1)))
            self.avelo1.append(self.avelo1[-1]-12*(math.pi**2)*\
            (self.x[-1]*math.sin(self.theta1[-1]) - self.y[-1]*math.cos(self.theta1[-1]))*\
            (self.x[-1]*math.cos(self.theta1[-1]) + self.y[-1]*math.sin(self.theta1[-1]))*\
            self.dt/(self.r[-1])**5)
            self.theta1.append(self.theta1[-1]+self.avelo1[-1]*self.dt)
            self.avelo2.append(self.avelo2[-1]-12*(math.pi**2)*\
            (self.x[-1]*math.sin(self.theta2[-1]) - self.y[-1]*math.cos(self.theta2[-1]))*\
            (self.x[-1]*math.cos(self.theta2[-1]) + self.y[-1]*math.sin(self.theta2[-1]))*\
            self.dt/(self.r[-1])**5)
            self.theta2.append(self.theta2[-1]+self.avelo2[-1]*self.dt)
            if self.theta1[-1] < -math.pi:
                self.theta1[-1] = self.theta1[-1] + 2*math.pi
            if self.theta1[-1] > math.pi:
                self.theta1[-1] = self.theta1[-1] - 2*math.pi
            else:
                pass
            if self.theta2[-1] < -math.pi:
                self.theta2[-1] = self.theta2[-1] + 2*math.pi
            if self.theta2[-1] > math.pi:
                self.theta2[-1] = self.theta2[-1] - 2*math.pi
            else:
                pass
            self.x.append(self.x[-1] + self.v_x[-1]*self.dt)
            self.y.append(self.y[-1] + self.v_y[-1]*self.dt)
            self.t.append(self.t[-1] + self.dt)
            self.dthe.append(abs(self.theta2[-1] - self.theta1[-1]))
    def show_results(self):
        #font = {'family': 'serif',
               #'color':  'darkred',
                #'weight': 'normal',
                #'size': 16,
                #}
        #plot,= pl.plot(self.t, self.theta1)
        #plot, = pl.plot(self.t,self.avelo1)
        plot,= pl.plot(self.t, self.dthe)
        pl.semilogy(self.t,self.dthe)
        pl.xlim(0,10)
        #pl.ylabel('$\Theta$(radians)')
        #pl.ylabel('$\omega$(radians/yr)')
        pl.ylabel('$\Delta$ $\Theta$(radians)')
        pl.xlabel('time(yr)')
        #pl.title('Hyperion $\Theta$ versus time 096')
        #pl.title('Hyperion $\omega$ versus time 096')
        pl.title('Hyperion $\Delta$ $\Theta$ versus time 096')
        pl.legend([plot,],['Circular orbit',],loc = "best")
        #pl.legend([plot,],['Elliptical orbit',],loc = "best")
        #pl.text(self.t[600], self.theta1[600],'Circular orbit', fontdict = font) 
        pl.show()
a = Hyperion(0,0.01)
a.calculate()
a.show_results()

import math as m
import pylab as pl
#Determine the initial value 
def initial(a,e): 
    x0=a*(1+e) 
    y0=0 
    v_x0=0 
    v_y0=2*m.pi*m.sqrt((1-e)/(a*(1+e))) 
    return [x0,y0,v_x0,v_y0] 
def orbits(beta, e, theta_0): 
    i_M=initial(1, e) 
    x0=i_M[0] 
    x=[] 
    x.append(x0) 
    y0=i_M[1] 
    y=[] 
    y.append(y0) 
    v_x0=i_M[2] 
    v_x=[] 
    v_x.append(v_x0) 
    v_y0=i_M[3] 
    v_y=[] 
    v_y.append(v_y0) 
    r=[] 
    r.append(m.sqrt(x0**2+y0**2)) 
    time=8
    t = [0]
    dt=0.0001 
    #Hyperion motion
    theta = [theta_0]
    w =[0]
    for i in range(int(time/dt)): 
        v_x.append(v_x[i]-4*m.pi**2*x[i]/(r[i]**(beta+1))*dt) 
        x.append(x[i]+v_x[i+1]*dt) 
        v_y.append(v_y[i]-4*m.pi**2*y[i]/(r[i]**(beta+1))*dt) 
        y.append(y[i]+v_y[i+1]*dt) 
        r.append(m.sqrt(x[i+1]**2+y[i+1]**2)) 
        w.append(w[i] - (3 * 4 * m.pi ** 2 * (x[i] * m.sin(theta[i]) - y[i] * \
                        m.cos(theta[i])) * (x[i] * m.cos(theta[i]) + y[i] * \
                        m.sin(theta[i]))) * dt)
        theta.append(theta[i] + w[i + 1] * dt)
        if theta[-1] >= m.pi:
            theta[-1] = theta[-1] - 2 * m.pi
        if theta[-1] <= -m.pi:
            theta[-1] = theta[-1] + 2 * m.pi
        t.append(t[i] + dt)
    return [x,y,v_x,v_y,r,theta,w,t] 
def delta_theta(theta_1, theta_2):
    delta_theta = []
    for i in range(len(theta_1)):
        delta_theta.append(abs(theta_1[i] - theta_2[i]))
    return delta_theta
#The orbits of moon 
M5 = orbits(2.0, 0, 0.004)
theta_M5 = M5[5]
w_M5 = M5[6]
t_M5 = M5[7]
M6 = orbits(2.0, 0, 0.003)
theta_M6 = M6[5]
w_M6 = M6[6]
t_M6 = M6[7]
delta_theta_M = delta_theta(theta_M6, theta_M5)
everage_M = sum(delta_theta_M) / len(delta_theta_M)#Lyapunov exponent
print everage_M
Delta_theta = []
for i in range(len(delta_theta_M)):
    Delta_theta.append(m.exp(everage_M * t_M5[i]))
pl.figure(figsize=[10,5])
pl.subplot(121)
pl.title('Hyperion  $\Delta\\theta$  versus  time',fontsize=15)
pl.xlabel('time $(yr)$',fontsize=15)
pl.ylabel('$\Delta\\theta$ (radius)',fontsize=15)
pl.text(3, 0.06, 'Circular orbit', fontsize=15)
pl.plot(t_M5, delta_theta_M)
pl.plot(t_M5, Delta_theta)
pl.semilogy(0.0001, 0.1)
M7 = orbits(2.0, 0.28, 0.004)
theta_M7 = M7[5]
w_M7 = M7[6]
t_M7 = M7[7]
M8 = orbits(2.0, 0.28, 0.003)
theta_M8 = M8[5]
w_M8 = M8[6]
t_M8 = M8[7]
delta_theta_M = delta_theta(theta_M8, theta_M7)
everage_M = sum(delta_theta_M) / len(delta_theta_M)
Delta_theta = []
for i in range(len(delta_theta_M)):
    Delta_theta.append(m.exp(everage_M * t_M7[i]))
print everage_M
pl.subplot(122)
pl.title('Hyperion  $\Delta\\theta$  versus  time',fontsize=15)
pl.xlabel('time $(yr)$',fontsize=15)
pl.ylabel('$\Delta\\theta$ (radius)',fontsize=15)
pl.plot(t_M7, delta_theta_M)
pl.plot(t_M7, Delta_theta)
pl.semilogy(0.0001, 0.1)

import matplotlib.pyplot as plt
import numpy as np
import math
class hyperion:
    def __init__(self,GM=4*math.pi**2,dt=0.0001,time=10):
        self.GM=GM
        self.x=[[1],[1],[1],[1]]
        self.y=[[0],[0],[0],[0]]
        self.vx=[[0],[0],[0],[0]]
        self.vy=[[5.5],[5],[4],[3]]
        self.x1=[[1],[1],[1],[1]]
        self.y1=[[0],[0],[0],[0]]
        self.vx1=[[0],[0],[0],[0]]
        self.vy1=[[5.5],[5],[4],[3]]
        self.dt=dt
        self.time=time
        self.r=[[math.sqrt(self.x[0][0]**2+self.y[0][0]**2)],[math.sqrt(self.x[1][0]**2+self.y[1][0]**2)],[math.sqrt(self.x[2][0]**2+self.y[2][0]**2)],[math.sqrt(self.x[3][0]**2+self.y[3][0]**2)]]
        self.t=[[0],[0],[0],[0]]
        self.w=[[0],[0],[0],[0]]
        self.theta=[[0.05],[0.05],[0.05],[0.05]]
        self.r1=[[math.sqrt(self.x1[0][0]**2+self.y1[0][0]**2)],[math.sqrt(self.x1[1][0]**2+self.y1[1][0]**2)],[math.sqrt(self.x1[2][0]**2+self.y1[2][0]**2)],[math.sqrt(self.x1[3][0]**2+self.y1[3][0]**2)]]
        self.t1=[[0],[0],[0],[0]]
        self.w1=[[0],[0],[0],[0]]
        self.theta1=[[0],[0],[0],[0]]
        self.d=[[math.log(abs(self.theta[0][0]-self.theta1[0][0]))],[math.log(abs(self.theta[1][0]-self.theta1[1][0]))],[math.log(abs(self.theta[2][0]-self.theta1[2][0]))],[math.log(abs(self.theta[3][0]-self.theta1[3][0]))]]
        self.dw=[[abs(self.w[0][0]-self.w1[0][0])],[abs(self.w[1][0]-self.w1[1][0])],[abs(self.w[2][0]-self.w1[2][0])],[abs(self.w[3][0]-self.w1[3][0])]]
    def calculate(self):
        for n in range(4):
            for i in range(int(self.time//self.dt)):
                self.vx[n].append(self.vx[n][i]-self.GM*self.x[n][i]*self.dt/self.r[n][i]**3)
                self.vy[n].append(self.vy[n][i]-self.GM*self.y[n][i]*self.dt/self.r[n][i]**3)
                self.x[n].append(self.x[n][i]+self.vx[n][i+1]*self.dt)
                self.y[n].append(self.y[n][i]+self.vy[n][i+1]*self.dt)
                self.r[n].append(math.sqrt(self.x[n][i+1]**2+self.y[n][i+1]**2))
                self.w[n].append(self.w[n][i]-3*self.GM/self.r[n][i]**5*(self.x[n][i]*math.sin(self.theta[n][i])-self.y[n][i]*math.cos(self.theta[n][i]))*(self.x[n][i]*math.cos(self.theta[n][i])+self.y[n][i]*math.sin(self.theta[n][i]))*self.dt)
                self.theta[n].append(self.theta[n][i]+self.w[n][i+1]*self.dt)
                self.t[n].append(self.t[n][i]+self.dt)
                if self.theta[n][i+1]<-math.pi:
                    self.theta[n][i+1]=self.theta[n][i+1]+2*math.pi
                if self.theta[n][i+1]>math.pi:
                    self.theta[n][i+1]=self.theta[n][i+1]-2*math.pi
                self.vx1[n].append(self.vx1[n][i]-self.GM*self.x1[n][i]*self.dt/self.r1[n][i]**3)
                self.vy1[n].append(self.vy1[n][i]-self.GM*self.y1[n][i]*self.dt/self.r1[n][i]**3)
                self.x1[n].append(self.x1[n][i]+self.vx1[n][i+1]*self.dt)
                self.y1[n].append(self.y1[n][i]+self.vy1[n][i+1]*self.dt)
                self.r1[n].append(math.sqrt(self.x1[n][i+1]**2+self.y1[n][i+1]**2))
                self.w1[n].append(self.w1[n][i]-3*self.GM/self.r1[n][i]**5*(self.x1[n][i]*math.sin(self.theta1[n][i])-self.y1[n][i]*math.cos(self.theta1[n][i]))*(self.x1[n][i]*math.cos(self.theta1[n][i])+self.y1[n][i]*math.sin(self.theta1[n][i]))*self.dt)
                self.theta1[n].append(self.theta1[n][i]+self.w1[n][i+1]*self.dt)
                self.t1[n].append(self.t1[n][i]+self.dt)
                if self.theta1[n][i+1]<-math.pi:
                    self.theta1[n][i+1]=self.theta1[n][i+1]+2*math.pi
                if self.theta1[n][i+1]>math.pi:
                    self.theta1[n][i+1]=self.theta1[n][i+1]-2*math.pi
                self.d[n].append(math.log(abs(self.theta[n][i]-self.theta1[n][i])))
                self.dw[n].append(abs(self.w[n][i]-self.w1[n][i]))
    def show_results(self,color):
        ax1=plt.subplot(221)
        ax2=plt.subplot(222)
        ax3=plt.subplot(223)
        ax4=plt.subplot(224)
        plt.sca(ax1)
        plt.plot(self.t[0],self.d[0],'g',label=r'd$\Theta$=0.01,v=5.5 Elliptical orbit')
        #plt.plot(self.t[0],self.dw[0],'b',label=r'd$\Theta$=0.05,v=5.5 Elliptical orbit')
        plt.title(r'Hyperion  $\Delta{\Theta}$ versus time',fontsize=12)
        #plt.title(r'Hyperion  $\Delta{\omega}$ versus time',fontsize=12)
        plt.xlabel(u'time(yr)',fontsize=12)
        plt.ylabel(u'$\Delta{\Theta}$(radians)',fontsize=12)
        #plt.ylabel(u'$\Delta{\omega}$(radians/yr)',fontsize=12)
        plt.legend(fontsize=10,loc='best')
        plt.sca(ax2)
        plt.plot(self.t[1],self.d[1],'g',label=r'd$\Theta$=0.01,v=5 Elliptical orbit')
        #plt.plot(self.t[1],self.dw[1],'b',label=r'd$\Theta$=0.05,v=5 Elliptical orbit')
        plt.title(r'Hyperion  $\Delta{\Theta}$ versus time',fontsize=12)
        #plt.title(r'Hyperion  $\Delta{\omega}$ versus time',fontsize=12)
        plt.xlabel(u'time(yr)',fontsize=12)
        plt.ylabel(u'$\Delta{\Theta}$(radians)',fontsize=12)
        #plt.ylabel(u'$\Delta{\omega}$(radians/yr)',fontsize=12)
        plt.legend(fontsize=10,loc='best')
        plt.sca(ax3)
        plt.plot(self.t[2],self.d[2],'g',label=r'd$\Theta$=0.01,v=4 Elliptical orbit')
        #plt.plot(self.t[2],self.dw[2],'b',label=r'd$\Theta$=0.05,v=4 Elliptical orbit')
        plt.title(r'Hyperion  $\Delta{\Theta}$ versus time',fontsize=12)
        #plt.title(r'Hyperion  $\Delta{\omega}$ versus time',fontsize=12)
        plt.xlabel(u'time(yr)',fontsize=12)
        plt.ylabel(u'$\Delta{\Theta}$(radians)',fontsize=12)
        #plt.ylabel(u'$\Delta{\omega}$(radians/yr)',fontsize=12)
        plt.legend(fontsize=10,loc='best')
        plt.sca(ax4)
        plt.plot(self.t[3],self.d[3],'g',label=r'd$\Theta$=0.01,v=3 Elliptical orbit')
        #plt.plot(self.t[3],self.dw[3],'b',label=r'd$\Theta$=0.05,v=3 Elliptical orbit')
        plt.title(r'Hyperion  $\Delta{\Theta}$ versus time',fontsize=12)
        #plt.title(r'Hyperion  $\Delta{\omega}$ versus time',fontsize=12)
        plt.xlabel(u'time(yr)',fontsize=12)
        plt.ylabel(u'$\Delta{\Theta}$(radians)',fontsize=12)
        #plt.ylabel(u'$\Delta{\omega}$(radians/yr)',fontsize=12)
        plt.legend(fontsize=10,loc='best')
a=hyperion()
a.calculate()
a.show_results('b')
plt.show()