MATLAB TUTORIAL, part 2.3: Pendulum

Pendulum

 

Pendulum with moving support


The length of the pendulum is increasing by stretching of the wire due to the weight of the bob. The effective spring constant for a wire of rest length \( \ell_0 \) is

\[ k = ES/ \ell_0 , \]
where the elastic modulus (Young's modulus) for steel is about \( E \approx 2.0 \times 10^{11} \) Pa and S is the cross-sectional area. With pendulum 3 m long, the static increase in elongation is about \( \Delta \ell = 1.6 \) mm, which is clearly not negligible. There is also dynamic stretching of the wire from the apparent centrifugal and Coriolis forces acting on the bob during motion. We can evaluate this effect by adapting a spring-pendulum system analysis. WE go from rectangular to polar coordinates:
\[ x= \ell \,\sin \theta = z_0 \left( 1 + \xi \right) \sin \theta , \qquad z= z_0 - \ell\,\cos \theta = z_0 \left[ 1 - \left( 1 + \xi \right) \cos \theta \right] , , \]
where \( z_0 = \ell_0 + mg/k \) is the static pendulum length, \( \ell = z_0 \left( 1 + \xi \right) \) is the dynamic length, ξ is the fractional string extension, and θ is the deflection angle. the equations of motion for small deflections are
\begin{eqnarray*} \left( 1 + \xi \right)\ddot{\theta} + 2\dot{\theta}\,\dot{\xi} + \omega_p^2 \theta &=& 0 , \\ \ddot{\xi} + \omega_s^2 \xi - \dot{\theta} + \frac{1}[2}\,\omega_p^2 \theta^2 &=& 0, \end{equarray*}
where the overdots denote differentiation with respect to time t, \( \omega_p^2 = g/z_0 \) is the square of the pendulum frequency, and \( \omega_s^2 = k/m \) is the spring (string) frequency, squared.

To get a feeling for how rigid and massive the pendulum support must be, we model the support as mass M kept in place by a spring of constant K, as shown in the picture below.

p = Rectangle[{-Pi/6 - 1.2, 2}, {-5, -5}]
a = Graphics[{Gray, p}]
coil = ParametricPlot[{t + 1.2*Sin[3*t], 1.2*Cos[3*t] - 1.2}, {t, -Pi/6 , 5*Pi - Pi/6}, FrameTicks -> None, PlotStyle -> Black]
line = Graphics[{Thick, Line[{{16.3843645, -1.2}, {20, -1.2}}]}]
r = Graphics[{EdgeForm[{Thick, Blue}], FaceForm[], Rectangle[{20, -5}, {26, 2}]}]
back = RegionPlot[-5 < x < 28 || -7 < y < -5, {x, -5, 28}, {y, -7, -5}, PlotStyle -> LightOrange]
line2 = Graphics[{Thick, Dashed, Line[{{23, -1}, {23, -14}}]}]
line3 = Graphics[{Thick, Line[{{23, -1}, {29, -11.5}}]}]
disk = Graphics[{Pink, Disk[{29.5, -12}, 1]}]
text = Graphics[Text["\[Theta]" , {25.1, -9.5}]]
Show[a, coil, line, r, back, line2, line3, disk, text]
Pendulum with moving support.

The natural frequency of the support is

\[ \Omega = \left( K/M \right)^{1/2} . \] The coupled equations for the system are
\begin{eqnarray*} \ddot{\theta} + \omega_0^2 \theta + \ddot{x}/\ell &=& 0 , \\ \left( 1 + m/M \right) \ddot{x} + \left( m/M \right) \ell \,\ddot{\theta} + \Omega^2 x &=& 0, \end{eqnarray*}
where x is the displacement of the rigid support. The former equation says that the effect of sway is to impart an additional angular acceleration \( - \ddot{x}/\ell \) to the pendulum for small angles of oscillation (otherwise, we have to use \( \omega_0^2 \sin \theta \) instead of \( \omega_0^2 \theta \) ).

 

Spring Pendulum



function main
tstart=0;
tstop=20;
z0=0; z0prime=1; theta0=0.1; theta0prime=0.1;
initial_w=[z0, theta0,z0prime,theta0prime]
k=50 ; m=2;  g=9.81;
l=0.5+m*g/k;
[times,sols]=ode45(@(t,w) diffeq(t,w,k,m,l,g), [tstart,tstop],initial_w)

figure
comet(sols(:,1),sols(:,2))
title('Theta vs z')
xlabel('z'),ylabel('theta (rad)')

figure
comet((times),sols(:,1))
title('z vs time')
xlabel('time'),ylabel('z')

figure
comet((times),sols(:,2))
title('Theta vs time')
xlabel('time'),ylabel('theta (rad)')

end
function dwdt=diffeq(t,w,k,m,l,g)
u1=w(1); u2=w(3); v1=w(2); v2=w(4);
du1dt=u2;
du2dt=-k/m*u1-g/(2*l)*v1^2+v2^2;
dv1dt=v2;
dv2dt=1/(1+2*u1)*(-g/l*v1-g/l*u1*v1-2*u2*v2);
dwdt=[du1dt;dv1dt;du2dt;dv2dt];


function spring

%Initial Conditions
initial_w=[0;0.1;1.5;0.1];
m=.2; k=5; g=9.8; l=.5+(m*g/k);

%Solves differential equation
[time,sols]=ode45(@(t,w) diffeq(t,w,m,k,g,l),[0,20],initial_w) %column 1 of solutions is z(t), column 2 is theta(t), column 3 is z'(t), column 4 is theta'(t)
%Plots path
figure
comet(sols(:,1),sols(:,2)) %Plots z vs theta
title('Displacement vs Theta')
xlabel('Theta (rad)')
ylabel('Displacement')

figure
comet(time,sols(:,1)) %Plots z vs time
title('Displacement vs Time')
xlabel('Time (s)')
ylabel('Displacement')

figure
comet(time,sols(:,2)) %Plots theta vs time
title('Theta vs Time')
xlabel('Time (s)')
ylabel('Theta (rad)')

end
function dwdt = diffeq(t,w,m,k,g,l)
z=w(1); theta=w(2); vz=w(3); vtheta=w(4) %Assigns values from vector w
dzdt=vz; dthetadt=vtheta;

dvzdt=-((k/m)*z + (g/(2*l))*theta^2 - vtheta^2)
dvthetadt= -((g/l)*theta + (g/(2*l))*z*theta + 2*vz*vtheta)/(1+2*z)

dwdt=[dzdt;dthetadt;dvzdt;dvthetadt]; %Puts derivatives into matrix dwdt

end

 

Double Pendulum


Double pendulum.
Double compound pendulum.
Animation of double pendulum.

A double pendulum consists of one pendulum attached to another. Double pendula are an example of a simple physical system which can exhibit chaotic behavior with a strong sensitivity to initial conditions. Several variants of the double pendulum may be considered; the two limbs may be of equal or unequal lengths and masses, they may be simple pendulums/pendula or compound pendulums (also called complex pendulums) and the motion may be in three dimensions or restricted to the vertical plane.

Consider a double bob pendulum with masses m1 and m2 attached by rigid massless wires of lengths \( \ell_1 \) and \( \ell_2 . \) Further, let the angles the two wires make with the vertical be denoted \( \theta_1 \) and \( \theta_2 , \) as illustrated at the left. Let the origin of the Cartesian coordinate system is taken to be at the point of suspension of the first pendulum, with vertical axis pointed up. Finally, let gravity be given by g. Then the positions of the bobs are given by

\begin{eqnarray*} x_1 &=& \ell_1 \sin \theta_1 , \\ y_1 &=& -\ell_1 \cos \theta_1 , \\ x_2 &=& \ell_1 \sin \theta_1 + \ell_2 \sin \theta_2 , \\ y_2 &=& -\ell_1 \cos \theta_1 - \ell_2 \cos \theta_2 . \end{eqnarray*}
The potential energy of the system is then given by
\[ \Pi = g\, m_1 \left( \ell_1 + y_1 \right) + g\,m_2 \left( \ell_1 + \ell_2 + y_2 \right) = g\,m_1 \left( 1- \cos \theta_1 \right) + g\, m_2 \left( \ell_1 + \ell_2 - \ell_1 \cos \theta_1 - \ell_2 \cos \theta_2 \right) , \]
and the kinetic energy by
\begin{eqnarray*} \mbox{K} &=& \frac{1}{2}\, m_1 v_1^2 + \frac{1}{2}\, m_2 v_2^2 \\ &=& \frac{1}{2}\, m_1 \ell_1^2 \dot{\theta}_1^2 + \frac{1}{2}\, m_2 \left[ \ell_1^2 \dot{\theta}_1^2 + \ell_2^2 \dot{\theta}_2^2 +2\ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) \right] \end{eqnarray*}
because
\[ v_1^2 = \dot{x}_1^2 + \dot{y}_1^2 = \ell_1^2 \dot{\theta}_1^2 , \qquad v_2^2 = \dot{x}_2^2 + \dot{y}_2^2 = \left( \ell_1 \cos\theta_1 \,\dot{\theta}_1 + \ell_2 \cos \theta_2 \,\dot{\theta}_2 \right)^2 + \left( \ell_1 \sin\theta_1 \,\dot{\theta}_1 + \ell_2 \sin \theta_2 \,\dot{\theta}_2 \right)^2 . \]
The Lagrangian is then
\begin{eqnarray*} {\cal L} &=& \mbox{K} - \Pi \\ &=& \frac{1}{2}\, m_1 \ell_1^2 \dot{\theta}_1^2 + \frac{1}{2}\, m_2 \left[ \ell_1^2 \dot{\theta}_1^2 + \ell_2^2 \dot{\theta}_2^2 +2\ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) \right] + \left( m_1 + m_2 \right) g\ell_1 \cos \theta_1 + m_2 g \ell_2 \cos \theta_2 . \end{eqnarray*}
Therefore, for θ1 we have
\begin{eqnarray*} \frac{\partial {\cal L}}{\partial \dot{\theta}_1} &=& m_1 \ell_1^2 \dot{\theta}_1 + m_2 \ell_1^2 \dot{\theta}_1 + m_2 \ell_1 \ell_2 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) , \\ \frac{{\text d}}{{\text d}t} \left( \frac{\partial {\cal L}}{\partial \dot{\theta}_1} \right) &=& \left( m_1 + m_2 \right) \ell_1^2 \ddot{\theta}_1 + m_2 \ell_1 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \ell_2 \dot{\theta}_2 \sin \left( \theta_1 - \theta_2 \right) \left( \dot{\theta}_1 - \dot{\theta}_2 \right) , \\ \frac{\partial {\cal L}}{\partial \theta_1} &=& - \ell_1 g \left( m_1 + m_2 \right) \sin \theta_1 - m_2 \ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \sin \left( \theta_1 - \theta_2 \right) . \end{eqnarray*}
So the Euler-Lagrange differential equation \( \frac{{\text d}}{{\text d}t} \left( \frac{\partial {\cal L}}{\partial \dot{q}} \right) = \frac{\partial {\cal L}}{\partial q} \) for θ1 becomes
\[ \left( m_1 + m_2 \right) \ell_1^2 \ddot{\theta}_1 + m_2 \ell_1 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) + m_2 \ell_1 \ell_2 \dot{\theta}_2^2 \sin \left( \theta_1 - \theta_2 \right) + \ell_1 g \left( m_1 + m_2 \right) \sin \theta_1 =0 . \]
Dividing through by \( \ell_1 , \) this simplifies to
\[ \left( m_1 + m_2 \right) \ell_1 \ddot{\theta}_1 + m_2 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) + m_2 \ell_2 \dot{\theta}_2^2 \sin \left( \theta_1 - \theta_2 \right) + g \left( m_1 + m_2 \right) \sin \theta_1 =0 . \]
Similarly, for θ2, we have
\begin{eqnarray*} \frac{\partial {\cal L}}{\partial \dot{\theta}_2} &=& m_2 \ell_2^2 \dot{\theta}_2 + m_2 \ell_1 \ell_2 \dot{\theta}_1 \cos \left( \theta_1 - \theta_2 \right) , \\ \frac{{\text d}}{{\text d}t} \left( \frac{\partial {\cal L}}{\partial \dot{\theta}_2} \right) &=& m_2 \ell_2^2 \ddot{\theta}_2 + m_2 \ell_1 \ell_2 \ddot{\theta}_1 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \ell_2 \dot{\theta}_1 \sin \left( \theta_1 - \theta_2 \right) \left( \dot{\theta}_1 - \dot{\theta}_2 \right) , \\ \frac{\partial {\cal L}}{\partial \theta_2} &=& \ell_1 \ell_2 m_2 \dot{\theta}_1 \dot{\theta}_2 \sin \left( \theta_1 - \theta_2 \right) - \ell_2 m_2 g\,\sin \theta_2 , \end{eqnarray*}
which leads to
\[ m_2 \ell_2 \ddot{\theta}_2 + m_2 \ell_1 \ddot{\theta}_1 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \dot{\theta}_1^2 \sin \left( \theta_1 - \theta_2 \right) + m_2 g \, \sin \theta_2 =0 . \]
The coupled second-order ordinary differential equations
\begin{eqnarray*} \left( m_1 + m_2 \right) \ell_1 \ddot{\theta}_1 + m_2 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) + m_2 \ell_2 \dot{\theta}_2^2 \sin \left( \theta_1 - \theta_2 \right) + g \left( m_1 + m_2 \right) \sin \theta_1 &=& 0 , \\ m_2 \ell_2 \ddot{\theta}_2 + m_2 \ell_1 \ddot{\theta}_1 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \dot{\theta}_1^2 \sin \left( \theta_1 - \theta_2 \right) + m_2 g \, \sin \theta_2 &=& 0 \end{eqnarray*}

The equations of motion can also be written in the Hamiltonian formalism. Computing the generalized momenta gives

\begin{eqnarray*} p_1 &=& \frac{\partial {\cal L}}{\partial \dot{\theta}_1} = \left( m_1 + m_2 \right) \ell_1^2 \dot{\theta}_1 + m_2 \ell_1 \ell_2 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) , \\ p_2 &=& \frac{\partial {\cal L}}{\partial \dot{\theta}_2} = m_2 \ell_2^2 \dot{\theta}_2 + m_2 \ell_1 \ell_2 \dot{\theta}_1 \cos \left( \theta_1 - \theta_2 \right) . \end{eqnarray*}
The Hamiltonian becomes
\[ H = \dot{\theta}_1 p_1 + \dot{\theta}_2 p_2 - {\cal L} . \]
For animationof double pendulum equations, see:
https://www.myphysicslab.com/pendulum/double-pendulum-en.html

function DoublePendulum
  [l1,l2,~,~,~,s]=pendata;  % parameters
  ax1=[-s s -s s]; % vector for axes
  T = 25; % Process duration, s
  v = 3;  % Speed of animation, positive integer: 1, 2, 3,...
% Initial values
  theta1_1 = pi/2;
  theta1 = 0;
  theta2_1 = pi/2;
  theta2 = 0;
  y0 = [theta1_1, theta1, theta2_1, theta2];

  [t,y] = ode45(@Equations, [0, T], y0);
  x1 = l1*sin(y(:,1));
  y1 = -l1*cos(y(:,1));
  x2 = x1 + l2*sin(y(:,3));
  y2 = y1 - l2*cos(y(:,3));

  figure
  tic
%   Animation
  for i = 1:v:numel(t)
   plot(0, 0,'k.',x1(i),y1(i),'b.',x2(i),y2(i),'r.','markersize',40);
   axis(ax1)
   line([0 x1(i)], [0 y1(i)],'Linewidth',2);
   line([x1(i) x2(i)], [y1(i) y2(i)],'linewidth',2,'color','r');
   xlabel(['\it\bf X \rm        Animation time, s:  ',num2str(toc,3)],...
     'HorizontalAlignment', 'left');
   ylabel('\it\bf Y');
   title(['Double Pendulum       \it t \rm= ',num2str(t(i),3)],...
     'fontsize',13, 'HorizontalAlignment', 'left', 'Position',[-.4*s s 0])
   drawnow
  end

%   Graphs
  figure
  plot(x1,y1,x2,y2,'r')
  xlabel('\it X','fontSize',12);
  ylabel('\it Y','fontSize',12);
  title('Chaotic Motion of a Double Pendulum')

  figure
  plot(t,y(:,1),t,y(:,3),'r','linewidth',2)
  axis([0,T,min(y(:,3)),max(y(:,3))])
  legend('\theta_1','\theta_2')
  xlabel('Time, s','fontSize',12)
  ylabel('\theta','fontSize',18,'fontweight','bold')
  title('\theta_1(0) = \pi/2  and  \theta_2(0) = \pi/2','fontsize',13)
 end

function z = Equations(~, y)
  [l1,l2,m1,m2,g,~]=pendata;
  z(4,1)=0;
  a = (m1+m2)*l1 ;
  b = m2*l2*cos(y(1)-y(3)) ;
  c = m2*l1*cos(y(1)-y(3)) ;
  d = m2*l2 ;
  e = -m2*l2*y(4)* y(4)*sin(y(1)-y(3))-g*(m1+m2)*sin(y(1)) ;
  f = m2*l1*y(2)*y(2)*sin(y(1)-y(3))-m2*g*sin(y(3)) ;

  z(1) = y(2);
  z(3) = y(4) ;
  z(2) = (e*d-b*f)/(a*d-c*b) ;
  z(4) = (a*f-c*e)/(a*d-c*b) ;
end

function [l1,l2,m1,m2,g,s] = pendata
  l1 = .3;  % length of the upper part of the double pendulum, m
  l2 = .3;  % length of the lower part of the double pendulum, m
  m1 = .1;  % mass of the upper part of the double pendulum, kg
  m2 = .1;  % mass of the lower part of the double pendulum, kg
  g = 9.8;  % gravitational acceleration, m/s^2
  s = l1+l2;
end



Complete