Chaotic Pendulum 

Andy Gijbels
andy.gijbels@student.kuleuven.be
gijbelsandy@hotmail.com
www.agshome.tk 

Copyright © 2006  by Andy Gijbels
All rights reserved 

Description 

This Maple Worksheet solves the "Chaotic Pendulum" problem in a clear and understanding way. The results are visualised in several graphs and simulations. 

 

Introduction 

A chaotic pendulum is a two-dimensional dynamical system. It consists of  a number of n rods connected to one an other by pivots and the pendulums contain point masses at the end of the light rods. The chaotic pendulum is an example of a physical system that exhibits chaotic behavior and shows a sensitive dependence on initial conditions. 

 

The equations derived for the motion of the chaotic pendulum are based on Kinematics and Newton's Laws. Since energy is conserved in this physical workout, the motion of the chaotic pendulum will continue indefinitely.  

Solving this problem with Maple allows varying the number of  pendulums to 5  without to long calculation periods, which is 'impossible' by hand. Maple also gives the possibility to visualise the results in a simulation. 

 

Parameters can be set in the initialisation section! 

 

Initialization 

Packages 

We will use pointplot from plots to draw the rods, disk from plottools to draw the masses and Norm from VectorCalculus to calculate the length of the velocity vectors. 

>	with(plots):with(plottools):with(VectorCalculus):

 

 

Setup 

Number of pendula 

>

n:=2;

 

 

Real time duration simulation in seconds 

>	Time:=20;                                          # seconds

 

 

Unable or disable trailfunction 

>	trail:=proc(i) if i=unable then w:=1; elif i=disable then w:=0 else error("You have to choose between enable or disable, or else the trailfunction will automatically be disabled"); end if end proc:w:=0:

 

Warning, `w` is implicitly declared local to procedure `trail` 

Possible values: unable or disable
Import message:  Unabling the trailfunction can cause long calculation periods!! 

>	w:=trail(unable);

 

 

 

Physical workout 

Parameters 

Gravitational acceleration 

>	g:=9.81;                                           #unit: m/s^2

 

 

Mass  

Standard all mass are taken equal, but you can change them manually if you wish. 

>	M:=[seq(1,m=1..n)];                                #unit: kg

 

 

 

Length 

Standard all lengths are taken equal, but you can change them manually if you wish. 

>	L:=[seq(1,m=1..n)];                                #unit: m

 

 

 

Initial angles  

Between the pendulums and the vertical imaginary line through the beginning points of the pendula. 

The initial angles can also be set differently.  

>	ang:=[seq(((2-i)*Pi/(2*i)+Pi/6),i=1..n)];          #unit: rad

 

 

 

Kinematics 

Positionvectors 

In function of the angles. 

The hang-up point is defined as the origin, the y-axis points vertically up and the x-axis points to the right. 

>	angles:=[seq(alpha[i],i=1..n)];

 

 

>	p:=array(1..n);

 

 

>	p[1]:=[L[1]*sin(angles[1](t)),-L[1]*cos(angles[1](t))];

 

 

>	for i from 2 to n do  p[i]:=[p[i-1][1]+L[i]*sin(angles[i](t)),p[i-1][2]-L[i]*cos(angles[i](t))]: end do;

 

 

 

Velocities 

Which are the derivates of the positions to the time. 

>	v:=array(1..n);

 

 

>	for i from 1 to n do   v[i]:=[diff(p[i][1],t),diff(p[i][2],t)]: end do;

 

 

 

 

Accelerations 

Which are the derivates of the velocities to the time. 

These equations we will need later to find a set of differential equations in function of the angles. 

>	a:=array(1..n);

 

 

>	for i from 1 to n do   a[i]:=[diff(v[i][1],t),diff(v[i][2],t)]: end do;

 

 

 

 

Forces 

Explanation 

We can derive two force equations ( for x- and y-direction)  for each pendulum. S[i] with i>1 stands for contact forces between the pendulums and S1 stands for the contact force between the hang-up and the first pendulum. 

 

Equations 

>	V:=array(1..n);

 

 

>	for i from 1 to n-1 do V[i]:=[M[i]*a[i][1]=-S[i]*sin(angles[i](t))+S[i+1]*sin(angles[i+1](t)),M[i]*a[i][2]=+S[i]*cos(angles[i](t))-S[i+1]*cos(angles[i+1](t))-M[i]*g]; end do;

 

 

>	V[n]:=[M[n]*a[n][1]=-S[n]*sin(angles[n](t)),M[n]*a[n][2]=+S[n]*cos(angles[n](t))-M[n]*g];

 

 

Eliminationprocedure contactforces 

Explanation 

We can't calculate the contact forces between the pendulums, so we will eliminate these forces out of our equations so they can't cause any trouble in later calculations. The elimination procedure will reduce the total amount of equations from 2xn to n. This does not form a problem because we only need n equations to find n parameters: the angles. 

 

Eliminations 

>	A:=array(1..n);

 

 

>	A[n]:=rhs(isolate(expand(V[n][1]*cos(angles[n](t))),S[n]*cos(angles[n](t))*sin(angles[n](t))))=rhs(isolate(expand(V[n][2]*sin(angles[n](t))),S[n]*cos(angles[n](t))*sin(angles[n](t)))):

 

>

for i from 1 to n-1 do #PART 1: ELIMINATION CONTACTFORCE NEXT PENDULUM loc1:=rhs(isolate(expand(V[n-i+1][1]),S[n-i+1]*sin(angles[n+1-i](t)))): V[n-i][1]:=subs(S[n+1-i]*sin(angles[n+1-i](t))=loc1,V[n-i][1]); loc2:=rhs(isolate(expand(V[n-i+1][2]),S[n-i+1]*cos(angles[n+1-i](t)))): V[n-i][2]:=subs(S[n+1-i]*cos(angles[n+1-i](t))=loc2,V[n-i][2]); #PART 2: ELIMINATION CONTACTFORCE PREVIOUS PENDULUM A[n-i]:=rhs(isolate(expand(V[n-i][1]*cos(angles[n-i](t))),S[n-i]*cos(angles[n-i](t))*sin(angles[n-i](t))))=rhs(isolate(expand(V[n-i][2]*sin(angles[n-i](t))),S[n-i]*cos(angles[n-i](t))*sin(angles[n-i](t)))); end do:

 

 

Differential equations 

Initialisation 

T is the time interval in seconds between two points in time where we will evaluate the angles the pendulums make. N_steps is the number of times we will evaluate. So N_steps x T is the real time duration of the simulation. 

>	T:=0.04;N_steps:=Time/T;

 

 

 

>	print("The real duration of the simulation is equal to",convert(N_steps*T,string),"s");

 

 

 

Initial conditions 

The initial angle velocities are set to zero. 

>	bgvw:=ListTools[Flatten](zip((i,j)->[j(0)=i,D(j)(0)=0],ang,angles))[];

 

 

 

Equations 

The total set of n differential equations. Notice that they can't be solved separately, because the motion of a particular pendulum depends on the other pendulums due to the contact forces. 

>	seq(A[i],i=1..n);

 

 

Numerical Solution 

We will solve the differential equations N_steps times with a time interval T between each solution. 

>	vgl:= dsolve({seq(A[i],i=1..n),bgvw},numeric,output=array([seq(m*T,m=0..N_steps-1)]),maxfun=500000):

 

>	Motionparameters:=convert(vgl[1,1],listlist);

 

 

>	data:=convert(vgl[2,1],listlist):

 

 

Conservation of Energy 

One Pendulum (for example the first pendulum) 

Kinetic Energy 

Ek = m.v²/2.
v has to be the absolute speed which is the Norm of the vector sum of the relative velocity (relative to the vertical line through the beginning point of the pendulum) and the absolute velocity of the previous pendulum. 

>	Ek:=array(1..500):

 

>	for theta from 1 to N_steps do Vs:=<0,0>; Vr:=<-L[1]*data[theta][3]*cos(data[theta][2]),-L[1]*data[theta][3]*sin(data[theta][2])>; Ek[theta]:=M[1]/2*Norm(Vr)**2; end do:

 

 

 

Potential Energy 

Ep = m.g.h. Because  the hang-up point is defined as the origin and the y-axis points vertically up, the height can be negative, so the potential energy can be negative too! 

>	Ep:=array(1..500):

 

>	for theta from 1 to N_steps do Epa:=array(1.. n); Ep[theta]:=-M[1]*g*L[1]*cos(data[theta][2]); end do:

 

Total Energy 

The total energy is the sum of the kinetic and potential energy and has to be constant if there is conservation of energy 

>

print("Total Energy = Brown,Kinetic Energy = green,Potential Energy = red");display(pointplot([seq([ito*T,Ep[ito]+Ek[ito]],ito=1..N_steps)],connect=true,color=brown),pointplot([seq([ito*T,Ep[ito]],ito=1..N_steps)],connect=true,color=red),pointplot([seq([ito*T,Ek[ito]],ito=1..N_steps)],connect=true,color=green));

 

 

 

Conclusion 

If the system contains more than one pendulum, we cannot speak of conservation of energy for each pendulum separately because of the interaction between the pendulums. 

 

Total System 

Kinetic Energy 

Ek = m.v²/2.
v has to be the absolute speed which is the Norm of the vector sum of the relative velocity (relative the the vertical line through the beginning point of the pendulum) and the absolute velocity of the previous pendulum. 

>	Ek:=array(1..500):

 

>

for theta from 1 to N_steps do Eka:=array(1.. n):  Vs:=<0,0>; Vr:=<-L[1]*data[theta][3]*cos(data[theta][2]),-L[1]*data[theta][3]*sin(data[theta][2])>; Eka[1]:=M[1]/2*Norm(Vr)**2; for j from 2 to n do Vs:=Vr+Vs: Vr:=<-L[1]*data[theta][3+2*(j-1)]*cos(data[theta][2+2*(j-1)]),-L[1]*data[theta][3+2*(j-1)]*sin(data[theta][2+2*(j-1)])>; Eka[j]:=M[j]/2*(Norm(Vs+Vr))**2; end do: Ek[theta]:=sum(Eka[m],m=1..n); end do:

 

 

 

Potential Energy 

Ep = m.g.h. Because  the hang-up point is defined as the origin and the y-axis points vertically up the height can be negative, so the potential energy can be negative too! 

>	Ep:=array(1..500):

 

>	for theta from 1 to N_steps do Epa:=array(1.. n); for j from 1 to n do Epa[j]:=-M[j]*g*sum(L[lambda]*cos(data[theta][lambda*2])   ,lambda=1..j); end do: Ep[theta]:=sum(Epa[zeta],zeta=1..n); end do:

 

Total Energy 

The total energy is the sum of the kinetic and potential energy and has to be constant if there is conservation of energy. 

>

print("Total Energy = Brown,Kinetic Energy = green,Potential Energy = red");display(pointplot([seq([ito*T,Ep[ito]+Ek[ito]],ito=1..N_steps)],connect=true,color=brown),pointplot([seq([ito*T,Ep[ito]],ito=1..N_steps)],connect=true,color=red),pointplot([seq([ito*T,Ek[ito]],ito=1..N_steps)],connect=true,color=green));

 

 

 

Conclusion 

The total amount of energy is constant, so we can conclude that there is conservation of energy for the entire system. 

 

Simulations 

Without trail 

Implementation 

>

SIM:=[]:

 

>	pos:=array(1..n);

 

 

>

for i from 1 to N_steps do lin:=seq(pointplot(([[0,0],[L[1]*sin(data[i][2]),-L[1]*cos(data[i][2])]]+[[sum(L[m]*sin(data[i][m*2]),m=1..j),sum(-L[m]*cos(data[i][m*2]),m=1..j)],[sum(L[m+1]*sin(data[i][m*2+2]),m=1..j),sum(-L[m+1]*cos(data[i][m*2+2]),m=1..j)]]),connect=true,thickness=2),j=1..n-1); dis:=display(disk([L[1]*sin(data[i][2]),-L[1]*cos(data[i][2])],M[1]/10,color=black),seq(display(disk([L[1]*sin(data[i][2]),-L[1]*cos(data[i][2])]+[sum(L[m+1]*sin(data[i][m*2+2]),m=1..j),sum(-L[m+1]*cos(data[i][m*2+2]),m=1..j)],M[j+1]/10,color=black)),j=1..n-1)); tek:=display({pointplot([[0,0],[L[1]*sin(data[i][2]),-L[1]*cos(data[i][2])]],connect=true,thickness=2),lin,dis,display(disk([0,0],0.04, color=black), scaling=constrained)}); SIM:=[op(SIM),tek]: end do:

 

Visualisation  

>	display(SIM,insequence=true,axes=NONE,scaling=constrained);

 

Please download workseet to see animation

 

 

With trail 

Implementation 

>	SIM:=[]:SPOOR:=[]:

 

>	pos:=array(1..n);

 

 

>

if w=1 then for i from 1 to N_steps do lin:=seq(pointplot(([[0,0],[L[1]*sin(data[i][2]),-L[1]*cos(data[i][2])]]+[[sum(L[m]*sin(data[i][m*2]),m=1..j),sum(-L[m]*cos(data[i][m*2]),m=1..j)],[sum(L[m+1]*sin(data[i][m*2+2]),m=1..j),sum(-L[m+1]*cos(data[i][m*2+2]),m=1..j)]]),connect=true,thickness=2),j=1..n-1); tek:=display({pointplot([[0,0],[L[1]*sin(data[i][2]),-L[1]*cos(data[i][2])]],connect=true,thickness=2),lin}); sp:=display({pointplot(([seq([L[1]*sin(data[lambda][2]),-L[1]*cos(data[lambda][2])]+[sum(L[m+1]*sin(data[lambda][m*2+2]),m=1..n-1),sum(-L[m+1]*cos(data[lambda][m*2+2]),m=1..n-1)],lambda=1..i)]),connect=true,color=black),tek}); SIM:=[op(SIM),tek]: SPOOR:=[op(SPOOR),sp]: end do: end if:

 

Visualisation 

>	if w=0 then print("The trail function is disabled, you can enable it in the initialisation section"); else  display(SPOOR,insequence=true,axes=NONE,scaling=constrained); end if;

 

Please download workseet to see animation

>

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.