张梓桐计算物理第七次作业

1.Problem

3.12，3.13，3.14

2.Abstract

Having dicussed the Euler method in calculating projectile movement of a cannon shell,in the passage,we are going to talk about the Euler-Cromer method of the better calculation of oscillatory under damped and external force.

3.Background

Now we start to figure out the movement rule of a damped and forced pendulum.With the basic knowledge of Newtoon second law and the definition of $w$, we have the following equation:

$$\begin{eqnarray*} \frac{d^{2}\theta}{dt^{2}}=-\frac{g}{l}sin\theta-q\frac{d\theta}{dt}+F_{D}sin(\Omega_Dt) \end{eqnarray*}$$
where $F_{D}$ is the module of the external force with $SI$ unit,and $\Omega_{D}$ is the frequency of external force,also with $SI$ unit.
Here we introduce the variables $w$ to transfer the secdonary differential equation into firs-order differential equation:
$$\begin{eqnarray*} \frac{dw}{dt}=-\frac{g}{l}sin\theta-q\frac{d\theta}{dt}+F_{D}sin(\Omega_Dt) \end{eqnarray*}$$
$$\begin{eqnarray*} \frac{d\theta}{dt}=w \end{eqnarray*}$$
Now we rewrite the physics formula into precodes,notcing the Euler-Cromer method,rather than Euler method,which will lead to a divergence in the infinity.

• For each time step $i$,calculate $w$ and $\theta$ at time step $i+1$.

• $w_{i+1}=w_{i}-[\frac{g}{l}sin\theta_{i}+qw_{i}-F_{D}sin(\Omega_Dt)]\Delta t$

• $\theta_{i+1}=\theta_{i}+w_{i+1}\Delta t$

• If $\theta_{i+1}$ is out of the $[-\pi,+\pi]$,add or substract $2\pi$ to keep it in this range.

• $t_{i+1}=t_{i}+\Delta t$

• Repeat for the deisred number of time steps.

In problem 3.13 and 3.14,we are asked to calculate the and compare the behacior of two,nearly identical pendulums.And we are going to draw a '$\Delta \theta$ versus time' plot.

3.Main body

3.1 Solution to 3.13，3.14

import math
import pylab as pl
class oscillatory_motion():
    def __init__(self):
        intial_theta=0.2
        intial_w=0
        time_of_duration=150
        self.theta=[intial_theta]
        self.t=[0]
        self.w=[intial_w]
        self.dt=0.04
        self.q=0.5
        self.F_D=1.2
        self.w_D=0.6666667
        self.nsteps = int(time_of_duration // self.dt + 1)
        self.w2=[intial_w]
        self.theta2=[intial_theta+0.001]
        self.a=[abs(self.theta[0]-self.theta2[0])]
    def run(self):
        for i in range(self.nsteps):#set an arbitrary time limit of the plot
            tmp_w=self.w[i]-(math.sin(self.theta[i])+self.q*self.w[i]-self.F_D*math.sin(self.w_D*self.t[i]))*self.dt
            self.w.append(tmp_w)  
            tmp_theta=self.theta[i]+self.w[i+1]*self.dt                        
            self.theta.append(tmp_theta)
           # while self.theta[-1]>1*math.pi:
            #    self.theta[-1]=self.theta[-1]-2*math.pi
            #while  self.theta[-1]<-1*math.pi:
             #   self.theta[-1]=self.theta[-1]+2*math.pi 
            tmp_w2=self.w2[i]-(math.sin(self.theta2[i])+self.q*self.w2[i]-self.F_D*math.sin(self.w_D*self.t[i]))*self.dt
            self.w2.append(tmp_w2)  
            tmp_theta2=self.theta2[i]+self.w2[i+1]*self.dt                        
            self.theta2.append(tmp_theta2)
           # while self.theta2[-1]>1*math.pi:
            #    self.theta2[-1]=self.theta[-1]-2*math.pi
            #while  self.theta2[-1]<-1*math.pi:
             #   self.theta2[-1]=self.theta2[-1]+2*math.pi                     
            self.t.append(self.t[i]+self.dt)
            self.a.append(abs(self.theta[i+1]-self.theta2[i+1]))
    def show(self):
        pl.semilogy(self.t, self.a,label='$F_{D}$=%.2f'%self.F_D)
        pl.xlabel('Time ($s$)')
        pl.ylabel('$\Delta \\theta $(radius)')
        pl.title('$\Delta \\theta $ versus time')
        pl.legend(loc='upper right',frameon = True)
        pl.grid(True)        
        pl.show()
a=oscillatory_motion()
a.run()
a.show()

3.2 Solution to 3.12

import math
import pylab as pl
class oscillatory_motion():
    def __init__(self):
        intial_theta=0.2
        intial_w=0
        time_of_duration=10000
        self.theta=[intial_theta]
        self.t=[0]
        self.w=[intial_w]
        self.dt=0.01
        self.q=0.5        
        self.F_D=1.2
        self.w_D=0.6666667
        self.nsteps = int(time_of_duration // self.dt + 1)
        self.tmp_w=[0]
        self.tmp_theta=[0]
    def run(self):
        for i in range(self.nsteps):#set an arbitrary time limit of the plot
            tmp_w=self.w[i]-(math.sin(self.theta[i])+self.q*self.w[i]-self.F_D*math.sin(self.w_D*self.t[i]))*self.dt
            self.w.append(tmp_w)  
            tmp_theta=self.theta[i]+self.w[i+1]*self.dt                        
            self.theta.append(tmp_theta)
            while self.theta[-1]>1*math.pi:
                self.theta[-1]=self.theta[-1]-2*math.pi
            while  self.theta[-1]<-1*math.pi:
                self.theta[-1]=self.theta[-1]+2*math.pi 
            self.t.append(self.t[i]+self.dt)
        for i in range(self.nsteps):           
           if self.t[i]%(2*math.pi/self.w_D)<0.01:
                self.tmp_w.append(self.w[i])
                self.tmp_theta.append(self.theta[i])
    def show(self):
        pl.plot(self.tmp_theta, self.tmp_w,'.',label='$F_{D}$=%.2f'%self.F_D)
        pl.xlabel('$\\theta$ (radians)')
        pl.ylabel('$w$(radius/s)')
        pl.title('$w$ versus $\\theta$')
        pl.legend(loc='upper right',frameon = True)
        pl.grid(True)        
        pl.show()
a=oscillatory_motion()
a.run()
a.show()

4.Result

4.1 Result with different intial angle.

Here the parameters are $g=l=9.8$, $\Omega_{D}=\frac{2}{3}$, $q=0.5$, $F_{D}=1.2$, $\Delta \theta =0.001$

Here are the expoential simulation plot:($\lambda=0.139$)

Here the parameters are $g=l=9.8$, $\Omega_{D}=\frac{2}{3}$, $q=0.5$, $F_{D}=0.5$, $\Delta \theta =0.001$

Here are the expoential simulation plot:($\lambda=-0.249$)

Here are the photos with different $F_{D}$ of '$\theta$ versus time' plots

Here are the photos with different $F_{D}$ of '$w$ versus time' plots

4.2 Result with different damping factors

Here the parameters are $g=l=9.8$, $\Omega_{D}=\frac{2}{3}$, $F_{D}=1.2$, $q=0.5$,$\Delta q =0.001$

Simulation：（$\lambda=0.189$）

Here the parameters are $g=l=9.8$, $\Omega_{D}=\frac{2}{3}$, $F_{D}=1.2$, $q=0.5$,$\Delta q =0.001$

Simulation plot（$\lambda=0$）：

4.3 Result with 3.12

Here the parameters are $g=l=9.8$, $\Omega_{D}=\frac{2}{3}$, $F_{D}=1.2$, $q=0.5$,$dt=0.001$,time duration =10000.This is a plot corresponding to a maximum of the frive force.

This is a plot corresponding to a $\frac{\pi}{4}$ phase of the frive force.

6.Ackanowledgement

After my dicsussion,here is one thing that I want to stress.Withoud the attention of the calculation of computer,I ignore the difference of fraction and decimal,which led to a serious and confusing problem,that is, I mistakely took $\frac{2}{3}$ as 0.67.This problem leads to a zero result with my external force and out of which I had struggled almost two hours to figure it out.And I want to thanks to Yukang,who found its out but he didn't know why at first and finally told by Wang shixing.Thanks to both of them.

7.Refernece

1.Yu kang for his help in $\frac{2}{3}$
2.Matplotilb
3.<< How to think like a computer scientist>>