第九次作业_计算程序

程序1

# basic packages needed
# matplotlib - for plot 2D lines
# numpy - for math
import matplotlib.pyplot as plt
import numpy as np
import math
# class CROMER
class CROMER_1(object):
    def __init__(self, _theta0=0.2, _omg0=0.0, _t0=0.0, _l=9.8,_g=9.8, _dt=0.04, _time=60.0,_q=0.5,_OMG_D=2.0/3):
        self.theta, self.omg, self.t = [_theta0], [_omg0], [_t0]
        self.l, self.g, self.dt, self.time, self.n= _l, _g, _dt, _time, int(_time/_dt)
        self.OMG_D=_OMG_D
        self.F_D=[0.0,0.5,1.2]
        self.I=range(self.n)
        self.q=[0,0.5,1.2]
    def calculate(self):
        for i in self.I:
            self.omg.append(self.omg[i]+(-self.g/self.l*math.sin(self.theta[i])-self.q[0]*self.omg[i]+self.F_D[1]*math.sin(self.OMG_D*self.t[i]))*self.dt)
            self.theta.append(self.theta[i]+self.omg[i+1]*self.dt)
            if self.theta[i+1]<-math.pi:
               self.theta[i+1]=self.theta[i+1]+2*math.pi
            if self.theta[i+1]>math.pi:
               self.theta[i+1]=self.theta[i+1]-2*math.pi
            self.t.append(self.t[i]+self.dt)
    def plot_theta(self,_ax):
        _ax.plot(self.t, self.theta, '--',label='q=0')
    def plot_w(self,_ax):
        _ax.plot(self.t,self.omg,'--',label='q=0')
class CROMER_2(CROMER_1):
    def calculate(self):
        for i in self.I:
            self.omg.append(self.omg[i]+(-self.g/self.l*math.sin(self.theta[i])-self.q[1]*self.omg[i]+self.F_D[1]*math.sin(self.OMG_D*self.t[i]))*self.dt)
            self.theta.append(self.theta[i]+self.omg[i+1]*self.dt)
            if self.theta[i+1]<-math.pi:
               self.theta[i+1]=self.theta[i+1]+2*math.pi
            if self.theta[i+1]>math.pi:
               self.theta[i+1]=self.theta[i+1]-2*math.pi
            self.t.append(self.t[i]+self.dt)
    def plot_theta(self,_ax):
        _ax.plot(self.t, self.theta, '--',label='q=0.5')
    def plot_w(self,_ax):
        _ax.plot(self.t,self.omg,'--',label='q=0.5')
class CROMER_3(CROMER_1):
    def calculate(self):
        for i in self.I:
            self.omg.append(self.omg[i]+(-self.g/self.l*math.sin(self.theta[i])-self.q[2]*self.omg[i]+self.F_D[1]*math.sin(self.OMG_D*self.t[i]))*self.dt)
            self.theta.append(self.theta[i]+self.omg[i+1]*self.dt)
            if self.theta[i+1]<-math.pi:
               self.theta[i+1]=self.theta[i+1]+2*math.pi
            if self.theta[i+1]>math.pi:
               self.theta[i+1]=self.theta[i+1]-2*math.pi
            self.t.append(self.t[i]+self.dt)
    def plot_theta(self,_ax):
        _ax.plot(self.t, self.theta, '--',label='q=1.2')
    def plot_w(self,_ax):
        _ax.plot(self.t,self.omg,'--',label='q=1.2')
# plot :
# EULER CROMER METHOD are used
fig= plt.figure(figsize=(14,7))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
comp= CROMER_1()
comp.calculate()
comp.plot_theta(ax1)
comp.plot_w(ax2)
comp=CROMER_2()
comp.calculate()
comp.plot_theta(ax1)
comp.plot_w(ax2)
comp=CROMER_3()
comp.calculate()
comp.plot_theta(ax1)
comp.plot_w(ax2)
ax1.set_title('Nonlinear Pendulum - Angle',fontsize=14)
ax2.set_title('Nonlinear Pendulum - W',fontsize=14)
ax1.set_xlabel('time(s)',fontsize=14)
ax1.set_ylabel('Angel (rad)',fontsize=14)
ax2.set_xlabel('time(s)',fontsize=14)
ax2.set_ylabel('w(rad/s)',fontsize=14)
ax1.legend(fontsize=12,loc='best')
ax2.legend(fontsize=12,loc='best')
plt.show(fig)

程序2

import matplotlib.pyplot as plt
import numpy as np
import math
theta0=0.2
theta1=0.2001
omg0=0.0
t0=0.0
l=9.8
g=9.8
dt=0.04
time=400.0
OMG_D=2.0/3
theta=[theta0]
theta1=[theta1]
omg1=[omg0]
omg=[omg0]
t1=[t0]
t=[t0]
t2=[t0]
n=int(time/dt)
OMG_D=OMG_D
F_D=[0.0,1.5,1.2]
I=range(n)
q=[0,1.5,0.5]
delta=[]
for i in I:
    omg.append(omg[i]+(-g/l*math.sin(theta[i])-q[2]*omg[i]+F_D[2]*math.sin(OMG_D*t[i]))*dt)
    theta.append(theta[i]+omg[i+1]*dt)
    if theta[i+1]<-math.pi:
       theta[i+1]=theta[i+1]+2*math.pi
    if theta[i+1]>math.pi:
       theta[i+1]=theta[i+1]-2*math.pi
    t.append(t[i]+dt)
for i in I:
    omg1.append(omg1[i]+(-g/l*math.sin(theta1[i])-q[2]*omg1[i]+F_D[2]*math.sin(OMG_D*t1[i]))*dt)
    theta1.append(theta1[i]+omg1[i+1]*dt)
    if theta1[i+1]<-math.pi:
       theta1[i+1]=theta1[i+1]+2*math.pi
    if theta1[i+1]>math.pi:
       theta1[i+1]=theta1[i+1]-2*math.pi
    t1.append(t1[i]+dt)
for i in I:
    l=theta1[i]-theta[i]
    #delta.append(l)
    if l>=0:
       l=l
       delta.append(l)
    if l<0:
       l=-l
    delta.append(l)
    t2.append(t2[i]+dt)
delta=delta[1:int(10000.00000)]
t1=t1[1:int(10000.00000)]
# EULER CROMER METHOD are used
plt.plot(t1, delta, '--',label='F_D=0.5')
plt.title(u'Nonlinear Pendulum - delta Angle',fontsize=14)
plt.xlabel(u'time(s)',fontsize=14)
plt.ylabel(u'delta angel (rad)',fontsize=14)
plt.legend(fontsize=12,loc='best')
plt.show()

程序3

# basic packages needed
# matplotlib - for plot 2D lines
# numpy - for math
import matplotlib.pyplot as plt
import numpy as np
import math
import pylab as pl
# class CROMER
class CROMER_1(object):
    def __init__(self, _theta0=0.2, _omg0=0.0, _t0=0.0, _l=9.8,_g=9.8, _dt=6.0*math.pi/1000, _time=6000.0,_q=0.5,_OMG_D=2.0/3):
        self.theta, self.omg, self.t = [_theta0], [_omg0], [_t0]
        self.l, self.g, self.dt, self.time, self.n= _l, _g, _dt, _time, int(_time/_dt)
        self.OMG_D=_OMG_D
        self.F_D=[0.0,0.5,1.2]
        self.I=range(self.n)
        self.q=[0.5,0.5,0.5]
        self.theta1=[]
        self.omg1=[]
        self.p=0
    def calculate(self):
        for i in self.I:
            self.omg.append(self.omg[i]+(-self.g/self.l*math.sin(self.theta[i])-self.q[0]*self.omg[i]+self.F_D[1]*math.sin(self.OMG_D*self.t[i]))*self.dt)
            self.theta.append(self.theta[i]+self.omg[i+1]*self.dt)
            if self.theta[i+1]<-math.pi:
               self.theta[i+1]=self.theta[i+1]+2*math.pi
            if self.theta[i+1]>math.pi:
               self.theta[i+1]=self.theta[i+1]-2*math.pi
            self.t.append(self.t[i]+self.dt)
            self.p=i%1000
            if self.p==0:
               self.theta1.append(self.theta[-1])
               self.omg1.append(self.omg[-1])
        self.theta1=self.theta1[100:]
        self.omg1=self.omg1[100:]
    def plot_theta(self,_ax):
        _ax.axes.scatter(self.theta1, self.omg1, marker='x',color='m',label='F_D=0.5')
        print len(self.theta1)
   # def plot_w(self,_ax):
       # _ax.plot(self.t,self.omg,'--',label='q=0')
class CROMER_2(CROMER_1):
    def calculate(self):
        for i in self.I:
            self.omg.append(self.omg[i]+(-self.g/self.l*math.sin(self.theta[i])-self.q[1]*self.omg[i]+self.F_D[2]*math.sin(self.OMG_D*self.t[i]))*self.dt)
            self.theta.append(self.theta[i]+self.omg[i+1]*self.dt)
            if self.theta[i+1]<-math.pi:
               self.theta[i+1]=self.theta[i+1]+2*math.pi
            if self.theta[i+1]>math.pi:
               self.theta[i+1]=self.theta[i+1]-2*math.pi
            self.t.append(self.t[i]+self.dt)
            if i%1000==0:
               self.theta1.append(self.theta[-1])
               self.omg1.append(self.omg[-1])
        self.theta1=self.theta1[100:]
        self.omg1=self.omg1[100:]
  # def plot_theta(self,_ax):
     #   _ax.plot(self.t, self.theta, '---',label='q=0.1')
    def plot_w(self,_ax):
        _ax.axes.scatter(self.theta1,self.omg1,marker='o',color='b',label='F_D=1.2')
class CROMER_3(CROMER_1):
    def calculate(self):
        for i in self.I:
            self.omg.append(self.omg[i]+(-self.g/self.l*math.sin(self.theta[i])-self.q[2]*self.omg[i]+self.F_D[1]*math.sin(self.OMG_D*self.t[i]))*self.dt)
            self.theta.append(self.theta[i]+self.omg[i+1]*self.dt)
            if self.theta[i+1]<-math.pi:
               self.theta[i+1]=self.theta[i+1]+2*math.pi
            if self.theta[i+1]>math.pi:
               self.theta[i+1]=self.theta[i+1]-2*math.pi
            self.t.append(self.t[i]+self.dt)
    def plot_theta(self,_ax):
        _ax.plot(self.t, self.theta, '--',label='q=0.5')
    def plot_w(self,_ax):
        _ax.plot(self.t,self.omg,'--',label='q=0.5')
# plot :
# EULER CROMER METHOD are used
fig= plt.figure(figsize=(14,7))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
comp= CROMER_1()
comp.calculate()
comp.plot_theta(ax1)
#comp.plot_w(ax2)
comp=CROMER_2()
comp.calculate()
#comp.plot_theta(ax1)
comp.plot_w(ax2)
#comp=CROMER_3()
#comp.calculate()
#comp.plot_theta(ax1)
#comp.plot_w(ax2)
ax1.set_title('Poincare section - $2n\pi$',fontsize=14)
ax2.set_title('Poincare section - $2n\pi$',fontsize=14)
ax1.set_xlabel('theta(rad)',fontsize=14)
ax1.set_ylabel('Angel volocity (rad/s)',fontsize=14)
ax2.set_xlabel('theta(rad)',fontsize=14)
ax2.set_ylabel('Angel volocity(rad/s)',fontsize=14)
ax1.legend(fontsize=12,loc='best')
ax2.legend(fontsize=12,loc='best')
plt.show(fig)