Consider a simplified model of the Newton’s cradle consisting of a pair of pendulums with mass as shown in Figure 1:

The Hybrid Model for the system is given by

where is the gravity constant, ; are the Pendulum’s moment of inertia and are the coefficients of viscous friction torque for each pendulum.

where is the mass ratio and is the restitution coefficient.

Source code for Simulation: script to run the simulation (run.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matlab M-file Project: HyEQ Toolbox @ Hybrid Dynamics and Control % Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software % % Filename: run % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function run %constants b_Ltilde = 0.25; b_Mtilde = 0.25; global gam l b_L b_R m_L m_R e; gam = 9.81; l = 1; m_L = 1; m_R = 2; b_L = b_Ltilde/m_L*l^2; b_R = b_Ltilde/m_R*l^2; e = 0.75; % initial conditions %x = [x1; x2; x3; x4] % x0 = [-pi/4; 0; pi/4; 0]; % case1 % x0 = [0;0; pi/4; 0]; % case2 x0 = [pi/4;0; pi/2; 0]; % case3 % simulation horizon TSPAN=[0 10]; JSPAN = [0 30]; % rule for jumps % rule = 1 -> priority for jumps % rule = 2 -> priority for flows rule = 1; options = odeset('RelTol',1e-4,'MaxStep',1e-3); % simulate [t j x] = HyEQsolver( @f,@g,@C,@D,... x0,TSPAN,JSPAN,rule,options); % state x1 = x(:,1); x2 = x(:,2); x3 = x(:,3); x4 = x(:,4); % plot solution figure(1) clf % hold on subplot(2,1,1),plotflows(t,j,x1) grid on ylabel('x1','rot',0) xlabel('t') axis([0 10 -1.75 1.75]) set(gca,'ytick',-pi/2:pi/4:pi/2) set(gca,'YTickLabel',{'-pi2','rem','0','rem','pi2'}) subplot(2,1,2),plotflows(t,j,x2) grid on ylabel('x2','rot',0) xlabel('t') axis([0 10 -3.8 3.8]) axis([0 10 -4.5 4.5]) %case 1 and 2 and 3 set(gca,'ytick',-pi:pi/2:pi) set(gca,'YTickLabel',{'-pi','rem','0','rem','pi'}) figure(2) clf % hold on subplot(2,1,1),plotflows(t,j,x3) grid on ylabel('x3','rot',0) xlabel('t') % axis([0 10 -1 1]) axis([0 10 -1.75 1.75]) set(gca,'ytick',-pi/2:pi/4:pi/2) set(gca,'YTickLabel',{'-pi2','rem','0','rem','pi2'}) subplot(2,1,2),plotflows(t,j,x4) grid on ylabel('x4','rot',0) xlabel('t') % axis([0 10 -2.5 2.5]) axis([0 10 -4.5 4.5]) %case 1 and 2 and 3 set(gca,'ytick',-pi:pi/2:pi) set(gca,'YTickLabel',{'-pi','rem','0','rem','pi'}) end % % plot hybrid arc % plotHybridArc(t,j,x) % xlabel('j') % ylabel('t') % zlabel('x')

Flow map (f.m)

function xdot = f(x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matlab M-file % Project: AME549_HW1_1.3 % % Name: f.m % % Description: Flow map % % Version: 1.0 % Required files: - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global gam l b_L b_R % state x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4); % flow map x1dot = x2; x2dot = -gam/l*sin(x1) - b_L*x2; x3dot = x4; x4dot = -gam/l*sin(x3) - b_R*x4; xdot = [x1dot; x2dot; x3dot; x4dot]; end

Flow set (C.m)

function [value discrete] = C(x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matlab M-file Project: HyEQ Toolbox @ Hybrid Dynamics and Control % Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software % % Filename: C.m % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Description: Flow set % Return C if outside of D, and 1 if inside C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % state x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4); if (x1 >= x3) && (x2 < x4) value = 1; elseif x1 < x3 value = 1; else value = 0; end end

Jump map (g.m)

function xplus = g(x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matlab M-file Project: HyEQ Toolbox @ Hybrid Dynamics and Control % Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software % % Filename: g.m % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Project: Simulation of a hybrid system (bouncing ball) % Description: Jump map %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global m_L m_R e % state x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4); % constant lam = m_L/(m_L + m_R); % jump map x1plus = x1; x2plus = (lam - (1-lam)*e)*x2 + (1-lam)*(1+e)*x4; x3plus = x3; x4plus = lam*(1+e)*x2 + (1-lam - lam*e)*x4; xplus = [x1plus; x2plus; x3plus; x4plus]; end

Jump set (D.m)

function xplus = g(x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matlab M-file Project: HyEQ Toolbox @ Hybrid Dynamics and Control % Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software % % Filename: g.m % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Project: Simulation of a hybrid system (bouncing ball) % Description: Jump map %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global m_L m_R e % state x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4); % constant lam = m_L/(m_L + m_R); % jump map x1plus = x1; x2plus = (lam - (1-lam)*e)*x2 + (1-lam)*(1+e)*x4; x3plus = x3; x4plus = lam*(1+e)*x2 + (1-lam - lam*e)*x4; xplus = [x1plus; x2plus; x3plus; x4plus]; end

Simulation Results are shown for and initial conditions:

Advertisements