NTRT Simulator  Version: Master
 All Classes Namespaces Files Functions Variables Typedefs Friends Pages
scaling_counterexample_nonlinear.m
1 % Scaling counterexample
2 % Drew Sabelhaus
3 % 7-13-14
4 
5 % This basic simulation is designed to show, by counterexample,
6 % why use of nondimensional numbers is necessary for NTRT.
7 % Equivalently, this would show why the commutative property of
8 % multiplication does not apply to scaling of units.
9 
10 clear all; clc; close all;
11 
12 % This disproves the following posulate:
13 % 10 * (length / time ^2) == (10 * length) / time^2
14 
15 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16 % A spring-mass-damper attached horizontally to a rod, hinged at the
17 % origin.
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19 p = {};
20 p.m = 1.5; % kg;
21 p.x_r = 1; % meters; NOTE THIS IS NOT USED FOR THE ANGULAR CASE.
22 p.l = 2; % meters;
23 p.g = 9.81; % m/s^2
24 p.k = 100; % N/m
25 p.c = 10; % N*s/m, or kg/sec
26 p.theta_0 = pi/4; % radians
27 p.theta_r = pi/3; % radians. This is for the rotational spring model.
28 
29 
30 % Forward-integrate (Euler style) for a handful of seconds
31 dt = 0.0001; % sec
32 t_final = 10; % sec
33 t = [0:dt:t_final];
34 num_timesteps = size(t,2);
35 y = zeros(2, num_timesteps); % because we have 2 state variables.
36 
37 % initialize. rod starts at angle theta_0, with zero velocity.
38 y(:,1) = [ p.theta_0; 0];
39 
40 for i=2:num_timesteps
41  % simulate.
42  theta = y(1,i-1);
43  theta_dot = y(2,i-1);
44  y(:, i) = ( rod_hinge_dynamics(theta, theta_dot, p)) * dt + y(:,i-1);
45 end
46 
47 hold on;
48 plot(t, y(1,:));
49 %plot(t, y(2,:), 'r');
50 axis([0 t_final 0.5 1.5]);
51 title('Original, unscaled simulation');
52 xlabel('time, sec');
53 ylabel('Radial position of rod, radians');
54 
55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 % Now, for comparison, scale gravity. This should show different results!
57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 
59 scaling_factor = 2;
60 g_scaled = scaling_factor * p.g;
61 p.g = g_scaled;
62 % We can use the same timestepping, but need to collect new data points:
63 y_new = zeros(2, num_timesteps);
64 y_new(:,1) = [ p.theta_0; 0];
65 
66 for i=2:num_timesteps
67  % simulate.
68  theta = y_new(1,i-1);
69  theta_dot = y_new(2,i-1);
70  y_new(:, i) = ( rod_hinge_dynamics(theta, theta_dot, p)) * dt + y_new(:,i-1);
71 end
72 
73 figure;
74 hold on;
75 plot(t, y_new(1,:));
76 %plot(t, y_new(2,:), 'r');
77 axis([0 t_final 0.5 1.5]);
78 title('With only gravity scaled');
79 xlabel('time, sec');
80 ylabel('Angular position of rod, radians');
81 
82 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83 % Here's the fun part: do a third simulation, and scale length,
84 % as we'd normally do in NTRT. Our length constants are the length of the
85 % rod, and the rest length of the spring.
86 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
87 
88 l_new = (1/scaling_factor) * p.l;
89 %l_new = scaling_factor * p.l;
90 %l_new = sqrt(scaling_factor) * p.l;
91 %x_r_new = scaling_factor * p.x_r;
92 %c_new = sqrt(scaling_factor) * p.c;
93 %c_new = (1/scaling_factor) * p.c;
94 p.l = l_new;
95 %p.x_r = x_r_new;
96 %p.c = c_new;
97 
98 %k_new = scaling_factor * p.k;
99 %m_new = (1/sqrt(scaling_factor)) * p.m;
100 %m_new = (1/scaling_factor) * p.m;
101 %p.k = k_new;
102 %p.m = m_new;
103 
104 % initialize.
105 y_new_length_scaled = zeros(2, num_timesteps);
106 y_new_length_scaled(:,1) = [p.theta_0; 0];
107 
108 for i=2:num_timesteps
109  % simulate.
110  theta = y_new_length_scaled(1,i-1);
111  theta_dot = y_new_length_scaled(2,i-1);
112  y_new_length_scaled(:, i) = ( rod_hinge_dynamics(theta, theta_dot, p)) * dt + y_new_length_scaled(:,i-1);
113 end
114 
115 
116 figure;
117 hold on;
118 plot(t, y_new_length_scaled(1,:));
119 %plot(t, y_new_length_scaled(2,:), 'r');
120 
121 %Adjust the limits of the plot to match our new scaling
122 %t_final_new = floor(t_final * (1/sqrt(scaling_factor)));
123 %axis([0 t_final 0.5 1.5]);
124 title('Gravity and length scale changed');
125 xlabel('time, sec');
126 ylabel('Angular position of rod, radians');
127 axis([0 t_final 0.5 1.5]);
128 
129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130 % Finally, scale according to our pi terms: both length and angular damping chang along with gravity.
131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132 
133 
134 %c_new = sqrt(scaling_factor) * p.c;
135 c_new = (1/scaling_factor) * p.c;
136 p.c = c_new;
137 
138 %k_new = scaling_factor * p.k;
139 %m_new = (1/sqrt(scaling_factor)) * p.m;
140 %m_new = (1/scaling_factor) * p.m;
141 %p.k = k_new;
142 %p.m = m_new;
143 
144 % initialize.
145 y_new_length_spring_scaled = zeros(2, num_timesteps);
146 y_new_length_spring_scaled(:,1) = [p.theta_0; 0];
147 
148 for i=2:num_timesteps
149  % simulate.
150  theta = y_new_length_spring_scaled(1,i-1);
151  theta_dot = y_new_length_spring_scaled(2,i-1);
152  y_new_length_spring_scaled(:, i) = ( rod_hinge_dynamics(theta, theta_dot, p)) * dt + y_new_length_spring_scaled(:,i-1);
153 end
154 
155 
156 figure;
157 hold on;
158 plot(t, y_new_length_spring_scaled(1,:));
159 title('Gravity, length, and angular damping constant scale changed');
160 xlabel('time, sec');
161 ylabel('Angular position of rod, radians');
162 axis([0 t_final 0.5 1.5]);
163 
164 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
165 % But wait! We forgot to do the time term.
166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
167 
168 
169 % initialize.
170 y_new_length_spring_time_scaled = zeros(2, num_timesteps);
171 y_new_length_spring_time_scaled(:,1) = [p.theta_0; 0];
172 
173 for i=2:num_timesteps
174  % simulate.
175  theta = y_new_length_spring_time_scaled(1,i-1);
176  theta_dot = y_new_length_spring_time_scaled(2,i-1);
177  y_new_length_spring_time_scaled(:, i) = ( rod_hinge_dynamics(theta, theta_dot, p)) * dt + y_new_length_spring_time_scaled(:,i-1);
178 end
179 
180 % Adjust time for the g, c, k angular scaling case
181 %t_new = t.* (1/scaling_factor);
182 t_new = t.* (scaling_factor);
183 
184 % Adjust the time scale, from our pi term analysis:
185 %t = t * scaling_factor;
186 
187 figure;
188 hold on;
189 plot(t_new, y_new_length_spring_time_scaled(1,:));
190 title('Gravity, length, angular damping constant, and time scale changed');
191 xlabel('time, sec');
192 ylabel('Angular position of rod, radians');
193 axis([0 t_final 0.5 1.5]);
194 
195 
void simulate(tgSimulation *simulation)