Rolling-Disk/Simulation.m at master · yardenas/Rolling-Disk · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
%% Lagrange and Newton derivation example comparison %%
omega0  = 3;   % [rad/sec]
R       = 0.5; % [m]
dtheta0 = 0.3; % [rad/sec] this is the small perturbation
IC = [0 0 0 omega0*R 0 0 0 dtheta0 0 omega0];

opt = odeset('reltol', 1e-8, 'abstol', 1e-8);
dt  = 1e-2;
[TimeL, XL] = ode45(@DerivativesL, 0:dt:15, IC, opt);
[TimeN, XN] = ode45(@DerivativesN, 0:dt:15, IC, opt);
plot(TimeL, XL(:,7),TimeN, XN(:,7),'--r');
xlabel('Time [sec]'); ylabel('\theta [rad]');
legend('Lagrange Derivation', 'Newton Derivation');
Animate(XL,dt);


%% Analytical stability derivation verification %%


%% 3Rdphi0^2 > g
clear all %#ok
close all
omega0  = 3;    % [rad/sec]
R       = 0.5;  % [m]
dtheta0 = 0.01; % [rad/sec] this is the small perturbation
dt  = 1e-2;
IC = [0 0 0 omega0*R 0 0 0 dtheta0 0 omega0];
opt = odeset('reltol', 1e-8, 'abstol', 1e-8);
[Time, X] = ode45(@DerivativesL, 0:dt:7, IC, opt);
figure();
plot(Time, X(:,7));
xlabel('Time [sec]'); ylabel('\theta [rad]');
Animate(X, dt);
figure()
th = X(:,7); dphi = X(:,end); dpsi = X(:,6);
m = 1; R = 0.5; alpha = 1/2;
for i = 1:length(X)
    dX(i,:) = DerivativesL(0,X(i,:).'); %#ok
end
ddphi = dX(:,end); ddpsi = dX(:,6); dth = X(:,8);
Const = sin(th).*(R*m*(R*ddphi + R*ddpsi.*sin(th) + 2*R*dpsi.*dth.*cos(th))...
    + R^2*alpha*m*(ddphi + ddpsi.*sin(th) + dpsi.*dth.*cos(th)))...
    + (R^2*alpha*m*cos(th).*(2*dphi.*dth + ddpsi.*cos(th)))/2;
plot(Time,Const);
xlabel('Time [sec]','interpreter','latex'); ylabel('$\dot{H_{a}}$','interpreter','latex');

%% 3Rdphi0^2 < g
clear all %#ok
close all
omega0  = 2;    % [rad/sec]
R       = 0.5;  % [m]
dtheta0 = 0.01; % [rad/sec] this is the small perturbation
IC = [0 0 0 omega0*R 0 0 0 dtheta0 0 omega0];
opt = odeset('reltol', 1e-8, 'abstol', 1e-8);
dt  = 1e-2;
[Time, X] = ode45(@DerivativesL, 0:dt:7, IC, opt);
figure
plot(Time, X(:,7));
xlabel('Time [sec]'); ylabel('\theta [rad]');
Animate(X, dt);


%% Analytical stability derivation verification (motion of a top) %%

%% 2.5Rdpsii0^2 > 2g
clear all %#ok
close all
g = 9.81; th = 10/180*pi; dpsi = 0.5; R = 0.5; psi = 0; alpha = 0.5;
dphi0 = -(4*g*sin(th) + 2*R*dpsi^2*sin(2*th)...
    + R*alpha*dpsi^2*sin(2*th))/(4*R*alpha*dpsi*cos(th)...
    + (4*R*dpsi*cos(psi)^2*cos(th))/(cos(psi)^2 + sin(psi)^2)...
    + (4*R*dpsi*cos(th)*sin(psi)^2)/(cos(psi)^2 + sin(psi)^2));

% IC = [0 0 0 0.5*dphi0 0 0.5 10/180*pi 0 0 dphi0]; %% circular patterns
IC = [0 0 0 0 0 6 0 0.1 0 0]; % stable top motion
dt = 1e-2;
opt = odeset('reltol', 1e-9, 'abstol', 1e-9);
[Time, X] = ode45(@DerivativesL, 0:dt:10, IC, opt);
figure()
plot(Time, X(:,7));
xlabel('Time [sec]'); ylabel('\theta [rad]');
Animate(X, dt);


%% 2.5Rdpsii0^2 < 2g
clear all %#ok
close all
IC = [0 0 0 0 0 3.5 0 0.1 0 0];
dt = 1e-2;
opt = odeset('reltol', 1e-9, 'abstol', 1e-9);
[Time, X] = ode45(@DerivativesL, 0:dt:10, IC, opt);
figure();
plot(Time, X(:,7));
xlabel('Time [sec]'); ylabel('\theta [rad]');
Animate(X,dt);