The 7th homework

Problem3.12 The chaos of pendulum

---

Abstract

This article mainly discusses chaos in thee driven nonlinear pendulum and involves the solutions of exercise 3.12.The crucial method used in the article is Euler-Cromer method.

---

background

●Problem3.12.
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≈2nπ/Ω, 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.

●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.

●In this paper, the relationship between angle and angular velocity is plotted, and the relationship between angular velocity and angle is drawn by simulating the method of flash display to calculate the chaotic phenomenon of nonlinear pendulum in physics textbooks. Figure, through these figures we can understand that chaotic phenomenon in the particle position is not completely random, but the average value with a certain law, specific to a certain point of relative fluctuations, these graphics can help us to compare the physical laws of intuitive Understanding, and thus appears to be crucial.

---

content

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 {g} {l}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 {g} {l}sinθ-qw+F_Dsin(Ω_Dt)$$
$$w=\frac {dθ} {dt}$$
he Euler-Cromer method is used for numerical calculation:
$$w_{i+1}=w_i-\frac {g} {l}sinθ△t-qw△t+F_Dsin(Ω_Dt)△t$$
$$θ_{i+1}=θ_i+w_i△t$$
$$t_{i+1}=t_i+△t$$
given properiate initial value, we can easily simulate the process. Here, we have shown several pictures of interest.

●The initial code is used to test the basic functions of the program

code1:

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()

result1:混沌现象1.png

●Instead of plotting θ as a function of t, let us plot the angular velocity w as a function of θ. It is called Phase-space plot.

code2

# -*- coding: utf-8 -*-
"""
@author: 胡胜勇
Last modified on  Mon Oct 31 00:10:24 2016
"""
# basic packages needed
#     matplotlib - for plot 2D lines
#     numpy - for math 
import matplotlib.pyplot as plt
import numpy as np
import math
# global var.
g = 9.8
l = 9.8
q = 0.5
OmgD = 2. / 3.
# class Euler
# use Euler METHOD to solve the pendulum
# para. & var.  :    
#             theta0, omg0, t0 - initial angel & angular velocity & time
#             l, g, dt, time - length of string, gravity acceleration, time step size & total time 
#             the time period of this pendulum - 1(s)
class Nonlinear_Cromer(object):
    """docstring for Nonlinear_Cromer"""
    global g, l, q, OmgD
    def __init__(self, _FD = .0, _theta0 = 0.2, _omg0 = .0, _t0 = .0, _dt = 0.04, _time = 400.):
        self.FD = _FD
        self.theta = [_theta0]
        self.omg = [_omg0]
        self.t = [_t0]
        self.dt = _dt
        self.time = _time
        self.n = int(_time / _dt)
        self.THETA = []
        self.OMG = []
        self.T =[]
# fun. calculate
# n=0 means reset
# n=1 means without reset   
    def calculate(self):
        for i in range(self.n):
            domg = (-(g / l) * math.sin(self.theta[-1]) - q * self.omg[-1] + self.FD * math.sin(OmgD * self.t[-1])) * self.dt
            self.omg.append(self.omg[-1] + domg)
            dtheta = self.omg[-1] * self.dt
            temp_theta = self.theta[-1] + dtheta
            if temp_theta < - math.pi:
                self.theta.append(temp_theta + 2 * math.pi)
            elif temp_theta > math.pi:
                self.theta.append(temp_theta - 2 * math.pi)
            else:
                self.theta.append(temp_theta)
            self.t.append(self.t[-1] + self.dt)
        # print 'theta length:',len(self.theta)
# fun. plotplot_phase_space
#  去除了阶跃点之间的连线的相空间图
    def plot_phase_space(self, _ax, _color):
        l = []
        for i, j in enumerate(self.theta):
            if abs(j - self.theta[i-1]) > 6.:
                l.append(i)
            else:
                pass
        l = [i -1 for i in l]
        l = [-1] + l + [len(self.theta)]
        # print l,len(l)
        for i, j in enumerate(l):
            if i < (len(l) - 1):
                self.THETA.append(self.theta[j+1 : l[i + 1]])
                self.OMG.append(self.omg[j+1 : l[i + 1]])
                self.T.append(self.t[j+1 : l[i + 1]])
        # print len(self.THETA), len(self.OMG), len(self.T)
        _ax.plot(self.THETA[0], self.OMG[0], '-', color = _color,label = r'$F_D = %.1f$' % self.FD)
        for i in range(len(self.T)-1):
            _ax.plot(self.THETA[i+1], self.OMG[i+1], '-', color = _color)
# plot:
#       ax1: theta vs time
fig = plt.figure(figsize = (16., 5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
comp = Nonlinear_Cromer(_FD = 0.5)
comp.calculate()
comp.plot_phase_space(ax1,'g')
comp = Nonlinear_Cromer(_FD = 1.2)
comp.calculate()
comp.plot_phase_space(ax2,'r')
for i in fig.axes:
    i.legend(fontsize=12,loc='upper right')
    # i.set_title(r'$\theta vs time$',fontsize=14)
    i.set_xlabel(r'$\theta (radians)$',fontsize=14)
    i.set_ylabel(r'$\omega (radians/s)$',fontsize=14)
    i.set_title(r'$Phase\quad Plot$'+r'$(q=\frac{1}{2}, l=g=9.8, \Omega_D=\frac{2}{3})$',fontsize=18)
# ax1.set_title(r'$\theta\quad vs\quad time$'+r'$(q=\frac{1}{2}, l=g=9.8, \Omega_D=\frac{2}{3})$',fontsize=18)
plt.show()

result2:混沌现象2.png

●we can also plot when $φ=0,φ=π/4,φ=π/2,φ=π$

code3:

# -*- coding: utf-8 -*-
"""
@author: 胡胜勇
Last modified on  Mon Oct 31 01:20:08 2016
"""
# basic packages needed
#     matplotlib - for plot 2D lines
#     numpy - for math 
import matplotlib.pyplot as plt
import numpy as np
import math
# global var.
g = 9.8
l = 9.8
q = 0.5
OmgD = 2. / 3.
# class Euler
# use Euler METHOD to solve the pendulum
# para. & var.  :    
#             theta0, omg0, t0 - initial angel & angular velocity & time
#             l, g, dt, time - length of string, gravity acceleration, time step size & total time 
#             the time period of this pendulum - 1(s)
class Nonlinear_Cromer(object):
    """docstring for Nonlinear_Cromer"""
    global g, l, q, OmgD
    def __init__(self, _FD = .0, _theta0 = 0.2, _omg0 = .0, _t0 = .0, _dt = 0.04, _time = 8000.):
        self.FD = _FD
        self.theta = [_theta0]
        self.omg = [_omg0]
        self.t = [_t0]
        self.dt = _dt
        self.time = _time
        self.n = int(_time / _dt)
        self.THETA = []
        self.OMG = []
        self.T =[]
# fun. calculate
# n=0 means reset
# n=1 means without reset   
    def calculate(self):
        for i in range(self.n):
            domg = (-(g / l) * math.sin(self.theta[-1]) - q * self.omg[-1] + self.FD * math.sin(OmgD * self.t[-1])) * self.dt
            self.omg.append(self.omg[-1] + domg)
            dtheta = self.omg[-1] * self.dt
            temp_theta = self.theta[-1] + dtheta
            if temp_theta < - math.pi:
                self.theta.append(temp_theta + 2 * math.pi)
            elif temp_theta > math.pi:
                self.theta.append(temp_theta - 2 * math.pi)
            else:
                self.theta.append(temp_theta)
            self.t.append(self.t[-1] + self.dt)
        # print 'theta length:',len(self.theta)
# fun. plotplot_phase_space
#  去除了阶跃点之间的连线的相空间图
    def plot_phase_space(self, _ax, _phi,_lab):
        t = np.array([(i * OmgD - _phi) / np.pi / 2 for i in self.t])
        tempt = np.array([round(i) for i in t])
        tempt2 = t - tempt
        a = []
        for i, j in enumerate(tempt2):
            if abs(j) < 0.001:
                a.append(i)
        self.THETA = [self.theta[i] for i in a]
        self.OMG = [self.omg[i] for i in a]
        # print t, tempt, tempt2,a
        # print a
        _ax.plot(self.THETA, self.OMG, 'o',markersize = 3.3,color = 'r',label = r'$F_D = %.1f, \phi=$' % self.FD + _lab )
        # for i in range(len(self.T)-1):
        #     _ax.plot(self.THETA[i+1], self.OMG[i+1], '-', color = _color)
# plot:
#       ax1: theta vs time
fig = plt.figure(figsize = (9.708, 6))
ax1 = fig.add_subplot(221, xlim = (-4, 4))
ax2 = fig.add_subplot(222, xlim = (-4, 4))
ax3 = fig.add_subplot(223, xlim = (-4, 4))
ax4 = fig.add_subplot(224, xlim = (-4, 4))
comp = Nonlinear_Cromer(_FD = 1.2)
comp.calculate()
comp.plot_phase_space(ax1,0,'0')
comp = Nonlinear_Cromer(_FD = 1.2)
comp.calculate()
comp.plot_phase_space(ax2,np.pi / 4,r'$\pi/4$')
comp = Nonlinear_Cromer(_FD = 1.2)
comp.calculate()
comp.plot_phase_space(ax3,np.pi / 2,r'$\pi/2$')
comp = Nonlinear_Cromer(_FD = 1.2)
comp.calculate()
comp.plot_phase_space(ax4, np.pi,r'$\pi$')
for i in fig.axes:
    i.legend(fontsize=12,loc='best')
    # i.set_title(r'$\theta vs time$',fontsize=14)
    i.set_xlabel(r'$\theta (radians)$',fontsize=14)
    i.set_ylabel(r'$\omega (radians/s)$',fontsize=14)
    # i.set_title(r'$Poincare Section$'+r'$(q=\frac{1}{2}, l=g=9.8, \Omega_D=\frac{2}{3})$',fontsize=18)
ax1.text(2,1.2,r'$Poincare Section$'+r'$(q=\frac{1}{2}, l=g=9.8, \Omega_D=\frac{2}{3})$'+'\n       ' +r'$\Omega_Dt=2n\pi +\phi$',fontsize=18)
plt.show()

result3:混沌现象3.png

Conclusion

It can be seen from the above calculation results that when the amplitude of the driving force is F_D = 0.5, the relationship between w and θ is relatively stable, that is, when taking a certain time, can determine the angular velocity at a certain angle, And the driving force is increased to 1.2, the relationship between the angle and the angular velocity is negative at a certain time, and the relationship between the angle and the angular velocity can not be uniquely determined, but it can be seen that these specific moments The relationship between the angle and angular velocity is always in accordance with certain laws of change, not completely random, this is the nature of the phenomenon of chaos.
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.

Acknowledgement

references:
1: Cai Hao. https://github.com/caihao/computational_physics_whu.
2.Brad Miller, David Ranum, Jeffrey Elkner, Peter Wentworth, Allen B. Downey, Chris Meyers, and Dario Mitchell. How to think like a computer scientist—— Learning with Python: Interactive Edition 2.0.
3.Nicholas J.Giordano,Hisao Naknishi.Computational Physics(Second Edition).
[4]. Wikipedia contributors, "Bifurcation diagram," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/w/index.php?title=Bifurcation_diagram&oldid=680690451.