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

The Hybrid Model for the system is given by

$f(x) =\left[ \begin{array}{c} x_2\\ -\frac{\gamma}{l} \cdot \sin(x_1)-b_L \cdot x_2\\ x_4\\ -\frac{\gamma}{l} \cdot \sin(x_3)-b_R \cdot x_4\\ \end{array} \right]$

where $\gamma$ is the gravity constant, $b_L= \frac{\widetilde{b_L}}{m_L l^2}, b_R= \frac{\widetilde{b_R}}{m_R l^2}$; $m_L l^2=m_R l^2$ are the Pendulum’s moment of inertia and $\widetilde{b_L}, \widetilde{b_R}$ are the coefficients of viscous friction torque for each pendulum.

$C:=\{x: x_1

$g(x) =\left[ \begin{array}{c} x_1\\ (\lambda -(1-\lambda)e) x_2+( 1- \lambda) (1+e) x_4 \\ x_3\\ \lambda (1+e) x_2 + (1-\lambda-\lambda e) x_4 \\ \end{array} \right]$

where $\lambda=\frac{m_L}{m_L+m_R}$ is the mass ratio and $e \in (0,1)$ is the restitution coefficient.

$D:=\{x: x_1=x_3, x_2>x_4\}$

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 $\widetilde{b_L} =\widetilde{b_R}= 0.25, \gamma : 9.81 [m/s^2], l=1 [m], m_L = 1; m_R = 2; e = 0.75$ and initial conditions: $x_0=\begin{bmatrix} \pi/4 & 0 & \pi/2 & 0 \\[0.3em] \end{bmatrix}^T$