A = 0.06;
k_l = 26400;
m = 483;
f = @(t,x,Om,l,k_s,d) [ x(2); ...
-(2*k_s*(x(1)-(A*sin(Om*t))))* ...
(sqrt((l-d)^2 + (x(1)-(A*sin(Om*t)))^2) - l)/ ...
(m*(sqrt((l-d)^2 + (x(1)-(A*sin(Om*t)))^2))) ];
Om_array = linspace(0,20,21);
l_array = linspace(0.2,1,21);
[om_array, L_array] = meshgrid(Om_array, l_array);
d = linspace(-0.005, -0.03, 10);
Response_amp = zeros([size(Om_array), numel(d)]);
T = 150;
x0 = [0,0];
for k=1:numel(d)
for i=1:numel(Om_array)
for j=1:numel(l_array)
Om = om_array(i,j);
l = L_array(i,j);
k_s = -(k_l*(l-d))/(4*d);
[t, x] = ode45(@(t,x) f(t,x,Om,l,k_s,d(k)),[100,T],x0);
Response_amp(i,j,k) = (max(x(:,1)) - min(x(:,1)))/2;
end
end
end
figure(1);
ax = axes();
view(3);
hold(ax);
view([30 33]);
grid on
for i=1:size(Response_amp,3)
mesh(om_array/(2*pi),L_array,Response_amp(:,:,i));
end
hold on
d = -1;
Response_amp = zeros([size(Om_array), numel(d)]);
T = 150;
x0 = [0,0];
for k=1:numel(d)
for i=1:numel(Om_array)
for j=1:numel(l_array)
Om = om_array(i,j);
l = L_array(i,j);
k_s = -(k_l*(l-d))/(4*d);
[t, x] = ode45(@(t,x) f(t,x,Om,l,k_s,d(k)),[100,T],x0);
Response_amp(i,j,k) = (max(x(:,1)) - min(x(:,1)))/2;
end
end
end
ax = axes();
view(3);
hold(ax);
view([30 33]);
grid on
for i=1:size(Response_amp,3)
mesh(om_array/(2*pi),L_array,Response_amp(:,:,i));
end
xlabel('Frequency (Hz)')
ylabel('Length of the spring (m)')
zlabel('Response Amplitude (m)')
set(gca,'FontSize',15)
hold off
