computationalphysics

exercise07:problem3.13

abstract

solution to problem3.13 of computational physics

key point:

1.plot the trajectory of two nearly identical pendulums
2.plot $\Delta \theta$ versus time
3.make a qualitative estimate of the corresponding Lyapunov exponent

background

physical equations

$$\frac{d\omega}{dt}=-\frac{g}{l}sin(\theta)-q\frac{d\theta}{dt}+F_Dsin(\Omega_Dt)$$
$$\frac{d\theta}{dt}=\omega$$

Euler_cromer calculation

$$\omega_{i+1}=\omega_i+[-\frac{g}{l}sin(\theta_i)-q\frac{d\theta_i}{dt}+F_Dsin(\Omega_Dt_i)]\Delta t$$
$$\theta_{i+1}=\theta_i+\omega_{i+1}\Delta t$$
if $\theta_{i+1}$is out of the range [-$\pi$,$\pi$],add or subtract $2\pi$ to keep it in this range.

parameters

$$l=g=9.8,q=1/2,\Omega_D=2/3,\Delta t=0.02$$
$$initial\space condition:\theta(0)=0.2,\omega(0)=0,,all\space in\space SI$$

main body

click here to check the main code

1.plot the trajectory of two nearly identical pendulums

a=chaotic_pendulums()
a.pendulum()
a.show_result()

here is the plot

the picture was uploaded to github,if failed loading,click here

2.plot $\Delta\theta$ versus $time$

b=chaotic_pendulums()
b.pendulum1()
b.show_result1()

here is the plot

the picture was uploaded to github,if failed loading,click here

3.make a qualitative estimate of the corresponding Lyapunov exponent

c=chaotic_pendulums()
c.cal_d_theta()
c.show_result2()
z=np.polyfit(a.t,a.d_theta,1)
p=np.poly1d(z)
print(z)
print(p)

the linear equation:

$$log(\Delta\theta)=y=0.08508 x - 7.159$$
thus,the Lyapunov exponent is:
$$\lambda=0.085>0$$

conclusion

while $F_D=1.2$,it is a chaotic pendulum,and the lyapunov exponent is about 0.085.

acknowledgement