Ex_08 - chaos in the nonlinear, damped, driven pendulum(ex3.18/21)

---

Question

Investigate the bifurcation diagrams found for the pendulum with other values of the drive frequency and damping parameter. Warning: this can easily become an ambitious project!
实际就是观察改变相应的外力频率和阻尼系数后bifurcation图会有什么变化，包括了3.18中的三个不同力的大小引起的变化。

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

接下来我选择重力加速度和摆长均为9.8，阻尼系数为0.5，外力频率为2/3，2/3+0.0001，2/3+0.0002，时间间隔为0.04.,制摆的角度与时间的关系图。

from __future__ import division
from math import *
from pylab import *
from numpy import *
# define a function that given amplitude of force, gives angle,avelo and t
def change_amp(F,f):
    # define some constants
    q = 0.5
    l = 9.8
    g = 9.8
    dt = pi/100
    theta0 = 0.2
    omega0 = 0
    angle = [theta0]
    avelo = [omega0]
    t = [0]
    an=[]
    # move the pendulum
    while t[-1] <= 1200*pi:
        avelo_new = avelo[-1] - ((g/l)*sin(angle[-1]) + q*avelo[-1] - F*sin(f*t[-1]))*dt
        avelo.append(avelo_new)
        angle_new = angle[-1] + avelo[-1]*dt
        while angle_new > pi:
            angle_new -= 2*pi
        while angle_new < -pi:
            angle_new += 2*pi
        angle.append(angle_new)
        if t[-1]>=900*pi:
            if t[-1]%(3*pi)<=dt:
                an.append(angle_new)
        t_new = t[-1] + dt
        t.append(t_new)
    return an
###
fdd_1=list(linspace(1.35,1.5,100))
fd_1=[]
th_1=[]
for i in fdd_1:
    points=change_amp(i,2/3)
    for j in points:
        fd_1.append(i)
        th_1.append(j)
###
fdd_2=list(linspace(1.35,1.5,100))
fd_2=[]
th_2=[]
for i in fdd_2:
    points=change_amp(i,2/3+0.00001)
    for j in points:
        fd_2.append(i)
        th_2.append(j)
###
fdd_3=list(linspace(1.35,1.5,100))
fd_3=[]
th_3=[]
for i in fdd_3:
    points=change_amp(i,2/3+0.00002)
    for j in points:
        fd_3.append(i)
        th_3.append(j)
###
subplot(1,3,1)
scatter(fd_1,th_1,s=10)
grid(True)
xlim(1.35,1.48)
ylim(0,3)
title('f=2/3')
xlabel('FD')
ylabel('theta(radians)')
subplot(1,3,2)
scatter(fd_2,th_2,s=10)
grid(True)
xlim(1.35,1.48)
ylim(0,3)
title('f=(2/3)+0.00001')
xlabel('FD')
ylabel('theta(radians)')
subplot(1,3,3)
scatter(fd_3,th_3,s=10)
grid(True)
xlim(1.35,1.48)
ylim(0,3)
title('f=(2/3)+0.00002')
xlabel('FD')
ylabel('theta(radians)')
show()

●当f增加时，图像的结构没有发生变化，只是图形整体下移，且上面的点有所增加。

Clink to cheek the reasult

from __future__ import division
from math import *
from pylab import *
from numpy import *
# define a function that given amplitude of force, gives angle,avelo and t
def change_amp(F,q):
    # define some constants
    l = 9.8
    g = 9.8
    dt = pi/100
    f=2/3
    theta0 = 0.2
    omega0 = 0
    angle = [theta0]
    avelo = [omega0]
    t = [0]
    an=[]
    # move the pendulum
    while t[-1] <= 1200*pi:
        avelo_new = avelo[-1] - ((g/l)*sin(angle[-1]) + q*avelo[-1] - F*sin(f*t[-1]))*dt
        avelo.append(avelo_new)
        angle_new = angle[-1] + avelo[-1]*dt
        while angle_new > pi:
            angle_new -= 2*pi
        while angle_new < -pi:
            angle_new += 2*pi
        angle.append(angle_new)
        if t[-1]>=900*pi:
            if t[-1]%(3*pi)<=dt:
                an.append(angle_new)
        t_new = t[-1] + dt
        t.append(t_new)
    return an
###
fdd_1=list(linspace(1.35,1.5,100))
fd_1=[]
th_1=[]
for i in fdd_1:
    points=change_amp(i,0.5)
    for j in points:
        fd_1.append(i)
        th_1.append(j)
###
fdd_2=list(linspace(1.35,1.5,100))
fd_2=[]
th_2=[]
for i in fdd_2:
    points=change_amp(i,0.51)
    for j in points:
        fd_2.append(i)
        th_2.append(j)
###
fdd_3=list(linspace(1.35,1.5,100))
fd_3=[]
th_3=[]
for i in fdd_3:
    points=change_amp(i,0.52)
    for j in points:
        fd_3.append(i)
        th_3.append(j)
###
subplot(1,3,1)
scatter(fd_1,th_1,s=10)
grid(True)
xlim(1.35,1.48)
ylim(0,3)
title('q=0.5')
xlabel('FD')
ylabel('theta(radians)')
subplot(1,3,2)
scatter(fd_2,th_2,s=10)
grid(True)
xlim(1.35,1.48)
ylim(0,3)
title('q=0.51')
xlabel('FD')
ylabel('theta(radians)')
subplot(1,3,3)
scatter(fd_3,th_3,s=10)
grid(True)
xlim(1.35,1.48)
ylim(0,3)
title('q=0.52')
xlabel('FD')
ylabel('theta(radians)')
show()

●其次改变阻尼系数，令q=0.5,0.51,0.52，结果如图。由图可知，当阻尼系数有微小的变化时，图像整体右移，亦即推迟了混沌现象的出现，如果持续增大阻尼系数，混沌现象会随之消失，回归到周期运动。

Clink to cheek the reasult
3.18题中改变力的大小，分别为1.4,1.44,1.465，用VPython可观察入下

from __future__ import division
from visual import *
from math import *
# add some constants
q=0.5
l=9.8
g=9.8
f=2.0/3
dt=0.04
theta0=0.2
omega0=0
# add three ceilings
ceil1=box(pos=(-20,5,0),size=(5,0.5,2),material=materials.bricks)
ceil2=box(pos=(0,5,0),size=(5,0.5,2),material=materials.bricks)
ceil3=box(pos=(20,5,0),size=(5,0.5,2),material=materials.bricks)
ball1=sphere(pos=(-20+l*sin(theta0),5-l*cos(theta0),0),radius=1,material=materials.emissive,color=color.red)
ball2=sphere(pos=(l*sin(theta0),5-l*cos(theta0),0),radius=1,material=materials.emissive,color=color.yellow)
ball3=sphere(pos=(20+l*sin(theta0),5-l*cos(theta0),0),radius=1,material=materials.emissive,color=color.green)
ball1.omega=omega0
ball2.omega=omega0
ball3.omega=omega0
ball1.theta=theta0
ball2.theta=theta0
ball3.theta=theta0
ball1.t=0
ball2.t=0
ball3.t=0
ball1.F=1.4
ball2.F=1.44
ball3.F=1.465
string1=curve(pos=[ceil1.pos,ball1.pos],color=color.red)
string2=curve(pos=[ceil2.pos,ball2.pos],color=color.yellow)
string3=curve(pos=[ceil3.pos,ball3.pos],color=color.green)
po1=arrow(pos=ball1.pos,axis=(10*ball1.omega*cos(ball1.theta),10*ball1.omega*sin(ball1.theta),0),color=color.red)
po2=arrow(pos=ball2.pos,axis=(10*ball2.omega*cos(ball2.theta),10*ball2.omega*sin(ball2.theta),0),color=color.yellow)
po3=arrow(pos=ball3.pos,axis=(10*ball3.omega*cos(ball3.theta),10*ball3.omega*sin(ball3.theta),0),color=color.green)
while True:
    rate(100)
    ball1.omega=ball1.omega-((g/l)*sin(ball1.theta)+q*ball1.omega-ball1.F*sin(f*ball1.t))*dt
    ball2.omega=ball2.omega-((g/l)*sin(ball2.theta)+q*ball2.omega-ball2.F*sin(f*ball2.t))*dt
    ball3.omega=ball3.omega-((g/l)*sin(ball3.theta)+q*ball3.omega-ball3.F*sin(f*ball3.t))*dt
    angle1_new=ball1.theta+ball1.omega*dt
    while angle1_new>pi:
        angle1_new-=2*pi
    while angle1_new<-pi:
        angle1_new+=2*pi
    angle2_new=ball2.theta+ball2.omega*dt
    while angle2_new>pi:
        angle2_new-=2*pi
    while angle2_new<-pi:
        angle2_new+=2*pi
    angle3_new=ball3.theta+ball3.omega*dt
    while angle3_new>pi:
        angle3_new-=2*pi
    while angle3_new<-pi:
        angle3_new+=2*pi
    ball1.theta=angle1_new
    ball2.theta=angle2_new
    ball3.theta=angle3_new
    ball1.pos=(-20+l*sin(ball1.theta),5-l*cos(ball1.theta),0)
    ball2.pos=(l*sin(ball2.theta),5-l*cos(ball2.theta),0)
    ball3.pos=(20+l*sin(ball3.theta),5-l*cos(ball3.theta),0)
    string1.pos=[ceil1.pos,ball1.pos]
    string2.pos=[ceil2.pos,ball2.pos]
    string3.pos=[ceil3.pos,ball3.pos]
    po1.pos=ball1.pos
    po2.pos=ball2.pos
    po3.pos=ball3.pos
    po1.axis=(10*ball1.omega*cos(ball1.theta),10*ball1.omega*sin(ball1.theta),0)
    po2.axis=(10*ball2.omega*cos(ball2.theta),10*ball2.omega*sin(ball2.theta),0)
    po3.axis=(10*ball3.omega*cos(ball3.theta),10*ball3.omega*sin(ball3.theta),0)
    ball1.t=ball1.t+dt
    ball2.t=ball2.t+dt
    ball3.t=ball3.t+dt

Conclusion

利用python来解决混沌问题是一个很好的选择，也对混沌有了更深刻的理解，发现步长以及初始条件对混沌有着很大的影响。同时经过分析Bifurcation diagram，知道了更多的知识，也进一步对python有了更深刻的理解。
利用VPython可以更直观的观察混沌摆的现象，更直接的研究其他因素的影响。但是3.21题改变阻尼系数和频率还不是很会做vpython展示，之后再做补充。