ANALYSIS OF A SYSTEM OF DOUBLE PENDULUM

SUBMISSION DATE: 07/04/2016

PHY 473A END SEMESTER REPORT ANALYSIS OF A SYSTEM OF DOUBLE PENDULUM Sudarshana Laha 151169 Shaswata Chowdhury 151168

1. Aim Of The Project: To solve the equation of motion of a double pendulum and display the motion of the system using animation which explicitly shows the behaviour of the system at different instants of time, given different initial conditions and hence study its chaotic behaviour and analyse the normal modes.

2. Assumptions: The assumptions taken into account to solve the problem are as follows: The bobs of the pendulum are considered to be point masses.  The strings attached to the bobs are massless and inextensible.

3. Deriving the equations of motion using Lagrangian of the system:

Fig 1. Schematic double pendulum setup As we can see from the schematic diagram we have two bobs of masses m1 and m2 respectively ,attached to strings of length l1 and l2 respectively. The angles each pendulum makes with the vertical are Ө1 and Ө2.

1

The expressions for kinetic energy and potential energy in terms of the polar coordinates are, given we have: 

x1 = l1 sin 𝜃1 , y1 = −𝑙1 cos 𝜃1 ,



2 1 Kinetic energy = 2 (𝑚1 + 𝑚2 )𝑙1 2 θ1̇ +



Potential energy = −(𝑚1 − 𝑚2 )𝑔𝑙1 cos 𝜃1 − 𝑚2 𝑙2 𝑔 cos 𝜃2

where 𝜃1̇ =

𝑑𝜃1 𝑑𝑡

and 𝜃2̇ =

𝑑𝜃2 𝑑𝑡

x2 = l1 sin 𝜃1 + l2 sin 𝜃2 y2 = −l1 cos 𝜃1 − l2 cos 𝜃2 1 2

2 𝑚2 𝑙2 2 θ2̇ + 𝑚2 𝑙1 𝑙2 𝜃1̇ 𝜃2̇ cos(𝜃1 − 𝜃2 )

are the angular velocities of the two bobs respectively.

We can formulate the Lagrangian, L= K.E – P.E as follows: 𝐿=

1 1 2 2 2 𝑚 𝑙 2 𝜃 ̇ + 𝑚2 [𝑙1 2 𝜃1̇ + 𝑙2 2 𝜃2̇ + 2𝑙1 𝑙2 𝜃1̇ 𝜃2̇ cos(𝜃1 − 𝜃2 )] + (𝑚1 + 𝑚2 )𝑔𝑙1 cos 𝜃1 + 𝑚2 𝑔𝑙2 cos 𝜃2 2 11 1 2

From the Euler Lagrange’s Equations we get the following equations of motion:    

𝑑𝜃1 𝑑𝑡 𝑑𝜃2 𝑑𝑡 𝑑𝜔1 𝑑𝑡 𝑑𝜔2 𝑑𝑡

= 𝜔1 = 𝜔2 = =

−𝑚2 𝑙1 𝜔1 2 sin(𝜃1 −𝜃2 ) cos(𝜃1 −𝜃2 )+𝑚2 𝑔 sin 𝜃2 cos(𝜃1 −𝜃2 )−𝑚2 𝑙2 𝜔2 2 sin(𝜃1 −𝜃2 )−(𝑚1 +𝑚2 )𝑔 sin 𝜃1 (𝑚1 +𝑚2 )𝑙1 −𝑚2 𝑙1 cos(𝜃1 −𝜃2 )2 −𝑚2 𝑙2 𝜔2 2 sin(𝜃1 −𝜃2 ) cos(𝜃1 −𝜃2 )+(𝑚1 +𝑚2 )(𝑔 sin 𝜃1 cos(𝜃1 −𝜃2 )+𝑙1 𝜔1 2 sin(𝜃1 −𝜃2 )−𝑔 sin 𝜃2 ) (𝑚1 +𝑚2 )𝑙2 −𝑚2 𝑙2 cos(𝜃1 −𝜃2 )2

4. Python Algorithm to solve the equations of motion: We know that of all the various methods that can be used to solve ordinary differential equations the RK4(Runge-Kutta) scheme provides the most accurate result as compared to Euler and RK2(Runge-Kutta) scheme. A second order method is one in which the error per time step is on the order of t3 . A fourth order method is one in which the error per time step is on the order of t5 . When the time step is less than one, the error decreases dramatically for the higher order method. So we work using the RK4 method with time step h or dt = 0.01.

2

To solve the ode of the form 𝒚̇ = 𝒇(𝒕, 𝒚) the RK4 uses a total of four tentative steps, k1-k4, to calculate a slope based on a weighted average where,



𝑘1 = 𝑓(𝑡𝑛 , 𝑦𝑛 )



𝑘2 = 𝑓 (𝑡𝑛 + 2 , 𝑦𝑛 + 2 𝑘1 )



𝑘3 = 𝑓 (𝑡𝑛 + 2 , 𝑦𝑛 + 2 𝑘2 )



𝑘4 = 𝑓(𝑡𝑛 + ℎ, 𝑦𝑛 + ℎ𝑘3 )

ℎ

ℎ

ℎ

ℎ

𝒉

And thus obtaining 𝒚𝒏+𝟏 =𝒚𝒏 + 𝟔 (𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒 )

5. Analysis of the results obtained from execution of code: We can check the validity of the program by plotting the energy of the system with time and since we run the code using dt=0.01 ,the fluctuation in energy is of the order of 10-10 .Then we do the same using dt=0.001. For initial conditions,𝜃1 = 5.7° , 𝜃2 = 0° , 𝜃1̇ = 0, 𝜃2̇ = 0 , 𝑚1 = 𝑚2 = 3𝑘𝑔 , 𝑙1 = 𝑙2 = 4 𝑚.

1.

dt=0.01

2.

dt=0.001

Fig 2. Plot of energy with time Here is enlisted 6 sets of conditions with time step=0.01 sec :   

𝑡 = 10𝑠𝑒𝑐, 𝑚1 = 𝑚2 = 3𝑘𝑔, 𝑙1 = 𝑙2 = 4𝑚, 𝜃1 = 10° , 𝜃2 = 0° , 𝜃1̇ = 𝜃2̇ = 0𝑠𝑒𝑐 −1 ………1 𝑡 = 10𝑠𝑒𝑐, 𝑚1 = 𝑚2 = 3𝑘𝑔, 𝑙1 = 𝑙2 = 4𝑚, 𝜃1 = 30° , 𝜃2 = 0° , 𝜃1̇ = 𝜃2̇ = 0𝑠𝑒𝑐 −1 ….……2 𝑡 = 10𝑠𝑒𝑐, 𝑚1 = 𝑚2 = 3𝑘𝑔, 𝑙1 = 𝑙2 = 4𝑚, 𝜃1 = 90° , 𝜃2 = 0° , 𝜃1̇ = 𝜃2̇ = 0𝑠𝑒𝑐 −1 ……….3

3

  

𝑡 = 30𝑠𝑒𝑐, 𝑚1 = 2𝑘𝑔, 𝑚2 = 1𝑘𝑔, 𝑙1 = 3𝑚, 𝑙2 = 4𝑚, 𝜃1 = 5.7° , 𝜃2 = 11.5° , 𝜃1̇ = 𝜃2̇ = 0𝑠𝑒𝑐 −1 ………4 ° ° ̇ 𝑡 = 50𝑠𝑒𝑐, 𝑚1 = 2𝑘𝑔, 𝑚2 = 1𝑘𝑔, 𝑙1 = 3𝑚, 𝑙2 = 4𝑚, 𝜃1 = 5.7 , 𝜃2 = −11.5 , 𝜃1 = 𝜃2̇ = 0𝑠𝑒𝑐 −1 ………5 𝑡 = 40𝑠𝑒𝑐, 𝑚1 = 3𝑘𝑔, 𝑚2 = 2𝑘𝑔, 𝑙1 = 4𝑚, 𝑙2 = 3𝑚, 𝜃1 = −90° , 𝜃2 = 0° , 𝜃1̇ = 2𝑠𝑒𝑐 −1 , 𝜃2̇ = 0𝑠𝑒𝑐 −1 .……….6

Phase Space plots: The phase space plots corresponding to 1,2,4,5 have a definite pattern and the pendula oscillates in small oscillation regime whereas 3,6 shows tendency towards chaotic behaviour.

1.

3.

2.

4.

4

5.

6.

Fig 3. Phase space plots of bob 1 and bob 2 Lissajous Figures: For small initial conditions corresponding to 1,2,4,5 the system behaves linearly. The motion is determined by simple sine and cosine curves. This is evident from the Lissajous curves, which plot with 𝜃1 respect to 𝜃2 . As the initial angles increase in magnitude, the pendulum is non-linear and is no more just a superposition of the sine and cosine curves. This seems to be a way of qualitatively observing chaos in the system as we can for 3,6.

1.

4.

2.

5.

3.

6.

Fig 4. Lissajous Figures 5

Variation of 𝜽𝟏 and 𝜽𝟐 with time: As we see in the Lissajous figures similarly a certain periodicity is observed in small oscillation regions in the variation of angles 𝜃1 and 𝜃2 . Corresponding to conditions 1, 2, 4, 5 we notice periodicity which vanishes for conditions 3,6.

2.

3.

5.

4.

6.

Fig 5. Variation of 𝜃1 and 𝜃2 with time 6

Animation: Corresponding to the 6 sets of conditions as mentioned above we have screenshots from the animations . 1 ,2 ,4 , 5 have more repetitive nature being in small oscillations domain compared to 3, 6 which show chaotic tendency.

1. At t=10 s

4.

At t=3.5 s

2. At t=4.8 s

5. At t=3.0 s

3. At t=2.9 s

6.

At t=1.9 s

Fig 6. Snapshots at different time for enlisted sets of conditions

Normal Modes Analysis: For this purpose we chose to study the case of small oscillations using the parameters: 1. dt = 0.01 s, total time =30 s , 𝑚1 = 𝑚2 = 1 𝑘𝑔, 𝑙1 = 𝑙2 = 1 𝑚, 𝜃1 = 5.72° , 𝜃2 = 5.72° , 𝜃1̇ = 0 𝑠𝑒𝑐 −1 , 𝜃2̇ = 0 𝑠𝑒𝑐 −1 (the in phase mode) 7

2. dt = 0.01 s, total time =30 s , 𝑚1 = 𝑚2 = 1 𝑘𝑔, 𝑙1 = 𝑙2 = 1 𝑚, 𝜃1 = −5.72° , 𝜃2 = 5.72° , 𝜃1̇ = 0 𝑠𝑒𝑐 −1 , 𝜃2̇ = 0 𝑠𝑒𝑐 −1 (the out of phase mode)

1.

2.

Fig 7. Real part of Fourier Transform of 𝜃1 and 𝜃2 . The peaks correspond to normal modes. The higher peak corresponds to the fourier transform of 𝜃2 . Using values from Fig. 7_1 the in phase value where 𝜃1 and 𝜃2 are in phase corresponds to frequency, f1=11.62 and from Fig. 7_2 where 𝜃1 and 𝜃2 are out of phase corresponds to frequency f2=27.78 . Using f1=11.58 ,the time period is 𝑇𝑃1 = 2𝜋 𝑇𝑃1

30 𝑓1

= 2.59 and thus normal mode frequency ,𝜔1 =

= 2.425 which is nearly close to the theoretical value as obtained in the code,i.e, 2.395.

Similarly, using f2=27.56 , the other 𝜔2 =5.772 which is close to theoretical value 5.7844.

6. In a nutshell:  We penned down the equations of motion by using Euler-Lagrange equation,decoupled the equations and formed 4 sets of linear ordinary differential equations.  We used RK-4 method to solve for 𝜃1 , 𝜃2 , 𝜃1̇ , 𝜃2.̇ numerically.

 We put the above values in the expression for energy and plotted it as function of time to validate our solution.

 We have formulated our code in a manner to make it flexible to a user who inputs his/her choice of bob masses, string lengths, initial angles, initial angular velocities, the time step and the total time for the simulation to run. 8

 We have plotted phase space plots, Lissajous figures , evolution of angles with time and also an animation showing the trajectory of two bobs to realize the nature of the system’s behaviour is very sensitive to initial conditions. Thus we studied both the periodic and chaotic nature of the double pendula system.

 We further analyzed the normal modes of oscillation at very small angles domain and obtained the normal frequency comparing it with our theoretically obtained values.

8. Acknowledgements: We are very grateful to Prof. M.K Verma for his valuable guidance and instilling in us the enthusiasm to solve such rich dynamical systems using the language Python. We would also like to thank the Teaching Assistants Abhishek Kumar and Anando Chatterjee for helping us out whenever we needed help.

9. Bibliography:    

http://www.physics.usyd.edu.au/~wheat/dpend_html/ https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/14IVPs/rk/complete.html http://www2.ph.ed.ac.uk/~egardi/MfP3-Dynamics/Dynamics_lecture_16.pdf https://www.google.co.in/search?q=lissajous+figures+of+double+pendulum&espv =2&biw=1366&bih=667&source=lnms&tbm=isch&sa=X&ved=0ahUKEwi0_-ycmnLAhWXjo4KHZX0D9sQ_AUIBygC&dpr=1

9