第八次作业

计算物理

---

正文

1.公式推导

$$\frac{d\omega}{dt}=-\frac{g}{l}sin\theta-q\frac{d\theta}{dt}+F_D(sin\Omega_D t)$$
$$\frac{d\theta}{dt}=\omega$$

2.代码

from math import *
from pylab import *
# define a function that given amplitude of force, gives angle,avelo and t
def change_amp(F):
    q = 0.5
    l = 9.8
    g = 9.8
    f = 2/3
    dt = pi/100
    theta0 = 0.2
    omega0 = 0
    angle = [theta0]
    avelo = [omega0]
    t = [0]
    an=[theta0]
    av=[omega0] 
    # move the pendulum
    while t[-1] < 3000*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)
        t_new = t[-1] + dt
        t.append(t_new)
        if t[-1]%(3*pi) <=dt:
            an.append(angle_new)
            av.append(avelo_new)
    return an,av
ang = change_amp(1.2)[0]
ave = change_amp(1.2)[1]
scatter(ang,)
title('angular velocity versus angle when F=1.2')
xlabel('theta(radians)')
ylabel('omega(radians/s)')
grid(True)
show()

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):
    # define some constants
    q = 0.5
    l = 9.8
    g = 9.8
    f = 2/3+0.00002
    dt = math.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
fd1=list(linspace(1.35,1.5,100))
fd=[]
th=[]
for i in fd1:
    points=change_amp(i)
    for j in points:
        fd.append(i)
        th.append(j)
scatter(fd,th,s=10)
grid(True)
xlim(1.35,1.5)
ylim(0,3)
title('Bifurcation diagram')
xlabel('FD')
ylabel('theta(radians)')
show()

"""
Created on Sat Nov 12 16:21:22 2016
@author: dw
"""
from __future__ import division
from visual import *
from math import *
q=0.5
l=9.8
g=9.8
f=2/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=0
ball2.F=0.5
ball3.F=1.2
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

3结果

1. $$当F_D取不同值时，\theta随时间变化$$

2. $$\theta随F_D的变化,当f变化很小时$$

3.用Vpython画出动图

2结论

参数的微小改变，对混沌现象的产生有着较大的影响

3感谢吴雨桥同学的代码提供