Exercise_7: Oscillatory and Motion and Chaos

physics python 3.13&3.14

Abstract

• Using Euler-Cromer method to solve the oscillatory motion.

• Also using class function,and take more factors into consideratioon.

Background

• We can use Euler method to solve some easy questions, but to solve the problem of the oscillatory motion and chaos with the conservation of energy, we cannot solve the problems easily by Euler method. we should use Euler-Cromer method to find the answer.

• we know the Euler-Cromer method. Now I will take the Euler-Cromer method into the oscillatory motion and chaos calculation:
  Newton’s second law for ideal pendulum with small angle:
  

  $$\frac{d^2\theta}{dx^2}=-\frac{l}{g}\theta$$

Write the second-order equations as two firest-order differential equations:

$$\frac{d\omega}{dt}=-\frac{g}{l}\omega$$

$$\frac{d\theta}{dt}=\omega$$

Finite difference form with Euler-Cromer method:

$$\omega_{i+1}=\omega_{i}-(\frac{g}{l})\theta_{i}\Delta{t}$$

$$\theta_{i+1}=\theta_{i}+\omega_{i+1}\Delta{t}$$

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

To make the question more realistic,we can take some factors into consideration.

1. We do not assume the small-angle approximation, and thus do not expand $sin\theta$ term in $\theta$

2. We include friction of the form

3. We add to our model a sinusoidal driving force $F_Dsin(\Omega_Dt)$

Putting all of these ingredients together, we have the equation of motion:

$$\frac{d^2\theta}{dt^2}=-\frac{g}{l}\theta-q\frac{d\theta}{dt}+F_Dsin(\Omega_Dt)$$
We can rewrite the two-order differential equation now as two first-order differential equations and obtian:
$$\frac{d\omega}{dt}=-\frac{g}{l}\theta-q\frac{d\theta}{dt}+F_Dsin(\Omega_Dt)$$

$$\frac{d\theta}{dt}=\omega$$
Finite difference form with Euler-Cromer method:
$$\omega_{i+1}=\omega_i-[(g/l)sin\theta_i-q\omega_i+F_Dsin(\Omega_Dt_i)]\Delta{t}$$

$$\theta_{i+1}=\theta_i+\omega_{i+1}\Delta{t}$$

$$t_{i+1}=t_i+\Delta{t}$$
If $\theta_{i+1}$ is out of the range $[-\pi,\pi]$ ,add or subtract $2\pi$ to keep it in this range.

Main body

• Problem 3.13 The initial values of θ differed by 0.001 rad.

Codes

import pylab as pl
import math
class oscillatory:
    def __init__(self,g=10,l=10,q=0.5,F_D=0.5,O_D=2/3,time_step=0.04,\
    total_time=50,initial_theta=0.2,initial_omega=0,initial_omega1=0,\
    initial_theta1=0.201,q1=0.5):
        self.g=g
        self.l=l
        self.q=q
        self.q1=q1
        self.F=F_D
        self.O=O_D
        self.t=[0]
        self.initial_theta=initial_theta
        self.initial_theta1=initial_theta1
        self.initial_omega=initial_omega
        self.dt=time_step
        self.time= total_time
        self.omega= [initial_omega]
        self.omega1=[initial_omega1]
        self.theta= [initial_theta]
        self.theta1= [initial_theta1]
        self.nsteps=int(total_time//time_step+1)
        self.D=[0]
        self.e=[0]    
    def run(self):
        for i in range(self.nsteps):
            tmpo=self.omega[i]+(-1*(self.g/self.l)*math.sin(self.theta[i])-\
            self.q*self.omega[i]+self.F*math.sin(self.O*self.t[i]))*self.dt           
            tmpt=self.theta[i]+tmpo*self.dt
            while(tmpt<(-1*math.pi) or tmpt>math.pi):
                if tmpt<(-1*math.pi):
                   tmpt+=2*math.pi
                if tmpt>math.pi:
                   tmpt-=2*math.pi
            self.omega.append(tmpo)
            self.theta.append(tmpt)
            tmpo1=self.omega1[i]+(-1*(self.g/self.l)*math.sin(self.theta1[i])-\
            self.q1*self.omega1[i]+self.F*math.sin(self.O*self.t[i]))*self.dt
            tmpt1=self.theta1[i]+tmpo1*self.dt
            while(tmpt1<(-1*math.pi) or tmpt1>math.pi):
                if tmpt1<(-1*math.pi):
                   tmpt1+=2*math.pi
                if tmpt1>math.pi:
                   tmpt1-=2*math.pi
            self.omega1.append(tmpo1)
            self.theta1.append(tmpt1)
            self.t.append(self.t[i]+self.dt)
            tmpD=abs(tmpt-tmpt1)
            self.D.append(tmpD)
            self.e.append(0.001*math.exp(-0.247*self.t[i]))
    def show_results(self):
        font = {'family': 'serif',
                'color':  'darkred',
                'weight': 'normal',
                'size': 16,}
        pl.semilogy(self.t,self.D)        
        pl.semilogy(self.t,self.e,'--')
        pl.title(r'$\Delta\theta$ versus time', fontdict = font)
        pl.xlabel('time(s)')
        pl.ylabel(r'$\Delta\theta$(radians)')
        pl.legend((['$F_D$=0.5']))
        pl.show()
a = oscillatory()
a.run()
a.show_results()

Reslut1:
1. F=0.5

2. F=1.2

• Problem3.14 the initial value of q differed by 0.001 rad.

Codes

import pylab as pl 
import math 
class oscillatory: 
    def __init__(self,g=9.8,l=9.8,q=0.5,F_D=1.2,O_D=2/3,time_step=0.04,\
    total_time=150,initial_theta=0.2,initial_omega=0,initial_omega1=0,\
    initial_theta1=0.2,q1=0.501): 
        self.g=g 
        self.l=l 
        self.q=q 
        self.q1=q1 
        self.F=F_D 
        self.O=O_D 
        self.t=[0] 
        self.initial_theta=initial_theta 
        self.initial_theta1=initial_theta1 
        self.initial_omega=initial_omega 
        self.dt=time_step 
        self.time= total_time 
        self.omega= [initial_omega] 
        self.omega1=[initial_omega1] 
        self.theta= [initial_theta] 
        self.theta1= [initial_theta1] 
        self.nsteps=int(total_time//time_step+1) 
        self.D=[0] 
        self.e=[0] 
    def run(self): 
        for i in range(self.nsteps): 
            tmpo=self.omega[i]+(-1*(self.g/self.l)*math.sin(self.theta[i])-\
            self.q*self.omega[i]+self.F*math.sin(self.O*self.t[i]))*self.dt 
            tmpt=self.theta[i]+tmpo*self.dt 
            while(tmpt<(-1*math.pi) or tmpt>math.pi): 
                if tmpt<(-1*math.pi): 
                   tmpt+=2*math.pi 
                if tmpt>math.pi: 
                   tmpt-=2*math.pi 
            self.omega.append(tmpo) 
            self.theta.append(tmpt) 
            tmpo1=self.omega1[i]+(-1*(self.g/self.l)*math.sin(self.theta1[i])-\
            self.q1*self.omega1[i]+self.F*math.sin(self.O*self.t[i]))*self.dt 
            tmpt1=self.theta1[i]+tmpo1*self.dt 
            while(tmpt1<(-1*math.pi) or tmpt1>math.pi): 
                if tmpt1<(-1*math.pi): 
                   tmpt1+=2*math.pi 
                if tmpt1>math.pi: 
                   tmpt1-=2*math.pi 
            self.omega1.append(tmpo1) 
            self.theta1.append(tmpt1) 
            self.t.append(self.t[i]+self.dt) 
            tmpD=abs(tmpt-tmpt1) 
            self.D.append(tmpD) 
            self.e.append(0.001*math.exp(0.25*self.t[i])) 
    def show_results(self): 
        font = {'family': 'serif', 
                'color':  'darkred', 
                'weight': 'normal', 
                'size': 16,} 
        pl.semilogy(self.t,self.D) 
        pl.semilogy(self.t,self.e,'--') 
        pl.title(r'$\Delta\theta$ versus time', fontdict = font) 
        pl.xlabel('time(s)') 
        pl.ylabel(r'$\Delta\theta$(radians)') 
        pl.legend((['$F_D$=1.2'])) 
        pl.show() 
a = oscillatory() 
a.run() 
a.show_results() 

Reslut2:
F=0.5

F=1.2

Conlusion

• We can use Euler-Cromer method to solve the oscillatory motion problem which cannot easily solved by Euler method.

Acknowledgment

• Thank rf from hongyi.