Double Spring Pendulum 

Andy Gijbels 

andy.gijbels@student.kuleuven.be 

gijbelsandy@hotmail.com 

www.agshome.tk 

Belgium 

Copyright ? 2007  by Andy Gijbels 

All rights reserved 

Description 

 

This Maple Worksheet solves the "Double Spring Pendulum" problem in a clear and understanding way. The results are visualised in a simulation. 

 

Introduction 

 

A double spring pendulum is a two-dimensional dynamical system. It consists of  a number of two springs connected to one an other by pivots. The pendulums contain point masses at the end of the weightless springs. The double spring 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 double spring  pendulum are based on Kinematics and Newton's Laws. Since energy is conserved in this physical workout, the motion of the chaotic spring pendulum will continue indefinitely.  

Solving this problem with Maple allows us to visualise the results in a nice simulation. 

 

Parameters can be set in the initialisation section! Beginvalues of motion are setted in the section "Beginvalues". 

 

Initialization 

Packages 

 

>	restart:with(plots):with(plottools):

 

 

Physics 

 

Parameters 

 

Universal parameters 

 

>	g:=9.81:                                                        # gravitational constant

 

Ball parameters 

 

>	m:=[4,4]:                                                       # mass balls

 

Spring parameters 

 

>	k:=[40,40]:                                                     # stiffness

 

>	l:=[2,2]:                                                       # length

 

Time parameters
 

>	T:=0.125:                                                       # step length

 

>	N_steps:= 150:                                                  # number of steps

 

Differential equations 

 

Method 

 

The equations are derived by using Newton's  law F = m.a 

F is the sum of all external forces working on the masses: spring forces and gravitational forces. 

The law is expressed in the two directions x ( horizontal ) and y  ( vertical ) 

 

Mass 1 

 

X direction 

 

>	vgl1:=k[1]*(l[1]-sqrt(x1(t)**2+y1(t)**2))*(x1(t)/sqrt(x1(t)**2+y1(t)**2))-k[2]*(l[2]-sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2)))*(x2(t)-x1(t))/sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2))= m[1]*diff(x1(t),t$2);

 

(4.2.2.1.1)

 

Y direction 

 

>	vgl2:=k[1]*(l[1]-sqrt(x1(t)**2+y1(t)**2))*(y1(t)/sqrt(x1(t)**2+y1(t)**2))-m[2]*g-k[2]*(l[2]-sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2)))*(y2(t)-y1(t))/sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2))= m[1]*diff(y1(t),t$2);

 

(4.2.2.2.1)

 

Mass 2 

 

X direction 

 

>	vgl3:=-k[2]*(l[2]-sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2)))*(x1(t)-x2(t))/sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2))= m[2]*diff(x2(t),t$2);

 

(4.2.3.1.1)

 

 

Y direction 

 

>	vgl4:=-m[2]*g-k[2]*(l[2]-sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2)))*(y1(t)-y2(t))/sqrt(((x1(t)-x2(t))**2+(y1(t)-y2(t))**2))= m[2]*diff(y2(t),t$2);

 

(4.2.3.2.1)

 

Initial values 

 

Remark 

 

The fixed point is the origin of the frame 

 

Mass 1 

 

Position 

 

>	P1:=y1(0)=-1,x1(0)=1;

 

(4.3.2.1.1)

 

Velocity 

 

>	V1:=D(y1)(0)=0,D(x1)(0)=0;

 

(4.3.2.2.1)

 

Mass 2 

 

Position 

 

>	P2:=y2(0)=-2,x2(0)=3;

 

(4.3.3.1.1)

 

Velocity 

 

>	V2:=D(y2)(0)=0,D(x2)(0)=0;

 

(4.3.3.2.1)

 

Collect values 

 

>	bgvw:=P1,V1,P2,V2;

 

(4.3.4.1)

 

Numeric solution 

 

>	sol:= dsolve({vgl1,vgl2,vgl3,vgl4,bgvw},numeric,output=array([seq(T*i,i=0..N_steps)]),maxfun=500000):

 

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

 

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

 

Graphics spring 

 

>	spring_help:=proc ()print("[x1,x2,y1,y2,number_of_coils,size_spring,time]") end proc;

 

(5.1)

 

>

veer:=proc(x1,y1,x2,y2,n,f,time) local kappa,lambda,theta,A,B,C,D,h,a,b,c,SIM,i,spring_static,g,k,X1,X2,Y1,Y2: kappa:=2*n; #_____________________________________________________ spring_static:=proc(x1,x2,y1,y2,n,h) NumericEventHandler(division_by_zero=proc(operator,operands,defVal): infinity end proc); lambda:=sqrt((x1-x2)**2+(y1-y2)**2)/kappa; theta:=arctan((y1-y2)/(x1-x2)); if (evalf(x1)>= evalf(x2) and evalf(y1)>= evalf(y2) and evalf(theta) >=0) then theta:=theta+Pi;end if: if (evalf(x1)>= evalf(x2) and evalf(y1)<=evalf(y2) and evalf(theta) <=0)then theta:=theta+Pi;end if: a[5]:=[x1,y1]; a[4]:=[x1,y1],[x1+(lambda/2)*cos(theta),y1+(lambda/2)*sin(theta)]; a[8]:=[x1+(lambda/2)*cos(theta),y1+(lambda/2)*sin(theta)],[x1+(lambda)*cos(theta)-h*sin(theta),y1+(lambda)*sin(theta)+h*cos(theta)]; a[7]:=[x2-(lambda)*cos(theta)-h*sin(theta),y2-(lambda)*sin(theta)+h*cos(theta)],[x2-(lambda/2)*cos(theta),y2-(lambda/2)*sin(theta)]; a[6]:=[x2-(lambda/2)*cos(theta),y2-(lambda/2)*sin(theta)],[x2,y2]; a[3]:=[x2,y2]; A:=seq(line([x1+(3+2*i)*lambda*cos(theta)-h*sin(theta),y1+(3+2*i)*lambda*sin(theta)+h*cos(theta)],[x1+2*(i+1)*lambda*cos(theta)+h*sin(theta),y1+2*(i+1)*lambda*sin(theta)-h*cos(theta)]),i=0..(n-2)/2-1);B:=line(a[4]),line(a[8]),line(a[6]),line(a[7]); C:=pointplot(a[3],thickness=4),pointplot(a[5],thickness=4); D:= seq(line([x1+(1+2*i)*lambda*cos(theta)-h*sin(theta),y1+(1+2*i)*lambda*sin(theta)+h*cos(theta)],  [x1+2*(i+1)*lambda*cos(theta)+h*sin(theta),y1+2*(i+1)*lambda*sin(theta)-h*cos(theta)]),i=0..(n-   2)/2-1); display(A,B,C,D,scaling=constrained); end proc: #_____________________________________________________ X1:=unapply(x1,x): X2:=unapply(x2,x): Y1:=unapply(y1,x): Y2:=unapply(y2,x): #_____________________________________________________ SIM:=[]; for i from 0 to time do   SIM:=[op(SIM),spring_static(X1(i),X2(i),Y1(i),Y2(i),kappa,f)]: end do: #_____________________________________________________ display(SIM,insequence=true,axes=none,scaling=constrained); end proc:

 

Simulating 

 

Implementation 

 

>

SIM:=[]:

 

>

for i from 1 to N_steps do tek:=display(veer(0,0,data[(i-1)+1][2],data[(i-1)+1][6],5,0.5,0),veer(data[(i-1)+1][2],data[(i-1)+1][6],data[(i-1)+1][4],data[(i-1)+1][8],5,0.5,0),disk([data[(i-1)+1][2],data[(i-1)+1][6]],0.2,color=black),disk([data[(i-1)+1][4],data[(i-1)+1][8]],0.2,color=black)): SIM:=[op(SIM),tek]: end do:

 

Visualisation 

 

>	display(SIM,insequence=true);

 

 

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.