@LiuYongJie 2016-10-30T15:43:01.000000Z 字数 7263 阅读 478

Ex_07 (3.12 chaos)

刘永杰 2014301020094

content

In constructing the Poincare section in Figure 3.9 we plotted points only at times that were in phase with the drive force; that is, at times $t≈2πn/Ω$ , where n is an integer. At these values of t the driving force passed through zero. However, we could just as easily have chosen to make the plot at times corresponding to a maximum of the drive force, or at times π/4 out-of-phase with this force, etc. Construct the Poincare sections for these cases and compare them with Figure 3.9.

Abstract

As we have investigatd the linear driven pendulum in the last section, this time, we will research on something more complicated as the nonliear system, or refered to chaorotic system.

Introduction

Chaos theory is the field of study in mathematics that studies the behavior of dynamical systems that are highly sensitive to initial conditions—a response popularly referred to as the butterfly effect.Small differences in initial conditions (such as those due to rounding errors in numerical computation) yield widely diverging outcomes for such dynamical systems, rendering long-term prediction impossible in general. This happens even though these systems are deterministic, meaning that their future behavior is fully determined by their initial conditions, with no random elements involved. In other words, the deterministic nature of these systems does not make them predictable. This behavior is known as deterministic chaos, or simply chaos.

Theory and results

Now that we have a numerical method that is suitable for various of the simple pendulum, and armed with some understanding of what might or might not happen when dissipation, an external force, and/or nonlinearity is present, we are ready to take on a slightly more complicated and also more interesting situation. that is, we add all three ingredients which were previously discussed separately. Put all these ingredients together, we have the equation of motion

$\frac {d^2θ} {dt^2}=- \frac gl sin(θ)-q\frac {dθ} {dt}+F_Dsin(Ω_Dt)$

we call this model for a nonlinear, dapmped, driven pendulum, the physical pendulum. It contains some very rich and interesting behaviors. However, as the equation can not be solved analytically, we have to turn to a numerical method. Under that case, the relevent differentlial equations are

$\frac {dw} {dt}=- \frac gl sin(θ)-q\frac {dθ} {dt}+F_Dsin(Ω_Dt)$

$\frac {dθ} {dt}=w$
given properiate initial value, we can easily simulate the process. Here, we have shown several pictures of interest.

code

最初编写的代码,用于测试程序的基本功能

import matplotlib.pyplot as plt
from math import *
print 'Exercise 3.12 test V1.0'
print 'Designed by roach'
Constant_g=9.8
Constant_l=9.8
Constant_q=0.5
Constant_Omegad=2.0/3.0
dt = 0.001
c_t = []
c_theta = []
c_omega = []
def Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega):
    c_theta.append(Initial_theta)
    c_omega.append(Initial_omega)
    c_t.append(0.0)
    print 'Initial_theta =',c_theta[0],'rad    ',
    print 'Initial_omega =',c_omega[0],'rad/s    ',
    print 'F_d =',F_d,'N    ',
    for i in range(100000):
        c_omega.append(c_omega[i]-(sin(c_theta[i])+Constant_q*c_omega[i]+F_d*sin(Constant_Omegad*c_t[i]))*dt)
        c_theta.append(c_theta[i] + c_omega[i+1]*dt)
        c_t.append(c_t[i]+dt)
        if c_theta[-1]>pi:
            c_theta[-1]=c_theta[-1]-2*pi
        if c_theta[-1]<-pi:
            c_theta[-1]=c_theta[-1]+2*pi
    print 'Total steps =',i,'    ',
    print 'dt =',dt
    return 0
Initial_theta=0.2
Initial_omega=0.0
F_dlist=[0,0.5,1.2]
for m in range(3):
    c_t = None
    c_theta = None
    c_omega = None
    c_t = []
    c_theta = []
    c_omega = []
    F_d=F_dlist[m]
    F_dstr=str(F_d)
    Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega)
    plt.plot(c_t,c_theta,label='F_D = '+F_dstr+' (N)')
plt.title('3.12 test')
plt.xlabel('t (s)')
plt.ylabel('theta (rad)')
plt.legend()
plt.show()

clink to view reasult

import matplotlib.pyplot as plt
from math import *
print 'Exercise 3.12 test V1.0'
print 'Designed by roach'
Constant_g=9.8
Constant_l=9.8
Constant_q=0.5
Constant_Omegad=2.0/3.0
dt = 0.01
c_t = []
c_theta = []
c_omega = []
def Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega):
    c_theta.append(Initial_theta)
    c_omega.append(Initial_omega)
    c_t.append(0.0)
    print 'Initial_theta =',c_theta[0],'rad    ',
    print 'Initial_omega =',c_omega[0],'rad/s    ',
    print 'F_d =',F_d,'Force unit    ',
    for i in range(100000):
        c_omega.append(c_omega[i]-(sin(c_theta[i])+Constant_q*c_omega[i]+F_d*sin(Constant_Omegad*c_t[i]))*dt)
        c_theta.append(c_theta[i] + c_omega[i+1]*dt)
        c_t.append(c_t[i]+dt)
        if c_theta[-1]>pi:
            c_theta[-1]=c_theta[-1]-2*pi
        if c_theta[-1]<-pi:
            c_theta[-1]=c_theta[-1]+2*pi
    print 'Total steps =',i,'    ',
    print 'dt =',dt
    return 0
Initial_theta=0.2
Initial_omega=0.0
F_dlist=[1.2,0.5,1.2]
for m in range(1):
    c_t = None
    c_theta = None
    c_omega = None
    c_t = []
    c_theta = []
    c_omega = []
    F_d=F_dlist[m]
    Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega)
    plt.scatter(c_theta,c_omega,color='blue',s=1)
plt.title('3.12 test')
plt.xlabel('theta (rad)')
plt.ylabel('omega (rad/s)')
plt.text(3,2,'F_D = 1.2(N)')
plt.show()

clink to view reasult
进一步选出同相的点

import matplotlib.pyplot as plt
from math import *
print 'Exercise 3.12 V1.0'
print 'Designed by roach'
Constant_g=9.8
Constant_l=9.8
Constant_q=0.5
Constant_Omegad=2.0/3.0
dt = 0.01
c_t = []
c_theta = []
c_omega = []
def Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega):
    c_theta.append(Initial_theta)
    c_omega.append(Initial_omega)
    c_t.append(0.0)
    print 'Initial_theta =',c_theta[0],'rad    ',
    print 'Initial_omega =',c_omega[0],'rad/s    ',
    print 'F_d =',F_d,'N    ',
    for i in range(100000):
        c_omega.append(c_omega[i]-(sin(c_theta[i])+Constant_q*c_omega[i]+F_d*sin(Constant_Omegad*c_t[i]))*dt)
        c_theta.append(c_theta[i] + c_omega[i+1]*dt)
        c_t.append(c_t[i]+dt)
        if c_theta[-1]>pi:
            c_theta[-1]=c_theta[-1]-2*pi
        if c_theta[-1]<-pi:
            c_theta[-1]=c_theta[-1]+2*pi
    print 'Total steps =',i,'    ',
    print 'dt =',dt
    return 0
Initial_theta=0.2
Initial_omega=0.0
F_dlist=[0,0.5,1.2]
for m in range(2,3):
    c_t = None
    c_theta = None
    c_omega = None
    c_t = []
    c_theta = []
    c_omega = []
    F_d=F_dlist[m]
    F_dstr=str(F_d)
    c_theta1=[]
    c_omega1=[]
    c_t1=[]
    Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega)
    for n1 in range(300):    
        for i in range(1,100000):
            if c_t[i-1] < 2*pi*n1/Constant_Omegad and c_t[i] > 2*pi*n1/Constant_Omegad:
                c_theta1.append(c_theta[i])
                c_omega1.append(c_omega[i])
    plt.scatter(c_theta1,c_omega1,label='theta in phase F_D = 1.2(N)')
plt.title('3.12 test')
plt.xlabel('theta (rad)')
plt.ylabel('omega (rad/s)')
plt.legend()
plt.show()

clink to view reasult
3.找出相位超前π/4的点(时间间隔 $dt=\frac{π}{300}s$ )

import matplotlib.pyplot as plt
from math import *
print 'Exercise 3.12 V1.0'
print 'Designed by roach'
Constant_g=9.8
Constant_l=9.8
Constant_q=0.5
Constant_Omegad=2.0/3.0
dt = pi/300
c_t = []
c_theta = []
c_omega = []
def Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega):
    c_theta.append(Initial_theta)
    c_omega.append(Initial_omega)
    c_t.append(0.0)
    print 'Initial_theta =',c_theta[0],'rad    ',
    print 'Initial_omega =',c_omega[0],'rad/s    ',
    print 'F_d =',F_d,'N    ',
    for i in range(3000000):
        c_omega.append(c_omega[i]-(sin(c_theta[i])+Constant_q*c_omega[i]+F_d*sin(Constant_Omegad*c_t[i]))*dt)
        c_theta.append(c_theta[i] + c_omega[i+1]*dt)
        c_t.append(c_t[i]+dt)
        if c_theta[-1]>pi:
            c_theta[-1]=c_theta[-1]-2*pi
        if c_theta[-1]<-pi:
            c_theta[-1]=c_theta[-1]+2*pi
    print 'Total steps =',i,'    ',
    print 'dt =',dt
    return 0
Initial_theta=0.2
Initial_omega=0.0
F_dlist=[0,0.5,1.2]
for m in range(2,3):
    c_t = None
    c_theta = None
    c_omega = None
    c_t = []
    c_theta = []
    c_omega = []
    F_d=F_dlist[m]
    F_dstr=str(F_d)
    c_theta1=[]
    c_omega1=[]
    c_t1=[]
    Calculate(Initial_theta,Initial_omega,c_t,c_theta,c_omega) 
    i=225
    n=0
    while i<3000000:
            c_theta1.append(c_theta[i])
            c_omega1.append(c_omega[i])
            c_t1.append(c_t[i])
            i=i+900
            n=n+1
    print 'Total',len(c_theta1)
    plt.scatter(c_theta1,c_omega1,label='theta pi/4 out of phase F_D = 1.2(N)')
plt.title('3.12 test')
plt.xlabel('theta (rad)')
plt.ylabel('omega (rad/s)')
plt.legend()
plt.show()

clink to view reasult

Conclusion

the real physical model is usually a very complicated system, chaos can happend in various system. In this article, we have only shown the chaorotic behaviors in a nonlinear, damping, driveing pendulum, more are still to be seen.