How do I simplify my code?
显示 更早的评论
clc, clear all
A= 0.06;
k_l = 26400; %Linear stiffness
m = 483; %Mass
l =0.5;
d =-0.005;
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
Om_array = linspace(0,20,20); %in rad/s-1
l_array = linspace(0,1,20);
[om_array, L_array] = meshgrid(Om_array, l_array);
Response_amp = zeros(size(Om_array));
T = 150;
x0 = [0,0];
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); %Spring stiffness
f = @(t,x) [ 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))) ];
[t, x] = ode45(f,[100,T],x0);
Response_amp(i,j) = (max(x(:,1)) - min(x(:,1)))/2;
% xval(i) = Om/(2*pi) ;
end
end
%% plot
figure(1);
ax = axes();
view(3);
hold(ax);
view([30 33]);
grid on
mesh(om_array/(2*pi),L_array,Response_amp) ;
xlabel('Frequency (Hz)')
ylabel('Length of the spring (m)')
zlabel('Response Amplitude (m)')
set(gca,'FontSize',15)
% set(gca,'xtick',[])
% set(gca,'ytick',[])
% set(gca,'ztick',[])
%%
%l = linspace(0,1,40);
%b = max(max(Response_amp));
hold on
d =-0.01;
Om_array = linspace(0,20,20); %in rad/s-1
l_array = linspace(0,1,20);
[om_array, L_array] = meshgrid(Om_array, l_array);
Response_amp = zeros(size(Om_array));
T = 150;
x0 = [0,0];
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); %Spring stiffness
f = @(t,x) [ 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))) ];
[t, x] = ode45(f,[100,T],x0);
Response_amp(i,j) = (max(x(:,1)) - min(x(:,1)))/2;
% xval(i) = Om/(2*pi) ;
end
end
%% plot
figure(1);
ax = axes();
view(3);
hold(ax);
view([30 33]);
grid on
mesh(om_array/(2*pi),L_array,Response_amp) ;
xlabel('Frequency (Hz)')
ylabel('Length of the spring (m)')
zlabel('Response Amplitude (m)')
set(gca,'FontSize',15)
mesh(om_array/(2*pi),L_array,Response_amp) ;
hold off
Hi, all. This is my code and this shows 2 graphs when d= -0.005 and d=-0.01. However, I wish to simplify my code as it seems to be long-winded. Also, I wish to variy d from -0.005 to -0.03. Thanks for reading.
采纳的回答
Ameer Hamza
2020-5-6
I am not sure if it can be sped up since you are alreay using ode45. But following simplfies the code and extend it to use for multiple values of 'd'
A = 0.06;
k_l = 26400; %Linear stiffness
m = 483; %Mass
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,20); %in rad/s-1
l_array = linspace(0,1,20);
[om_array, L_array] = meshgrid(Om_array, l_array);
d = linspace(-0.005, -0.01, 5);
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); %Spring stiffness
[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
%% plot
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
xlabel('Frequency (Hz)')
ylabel('Length of the spring (m)')
zlabel('Response Amplitude (m)')
set(gca,'FontSize',15)
% set(gca,'xtick',[])
% set(gca,'ytick',[])
% set(gca,'ztick',[])
4 个评论
Donghun Lee
2020-5-12
Ameer Hamza
2020-5-12
You can use subplot to plot each surface on separate axes to make the differences more clear. Although they still look quite similar, you can spot minor differences
figure;
for i=1:size(Response_amp,3)
subplot(2, 3, i)
mesh(om_array/(2*pi),L_array,Response_amp(:,:,i));
view([30 33]);
grid on
title(sprintf('d=%.4f', d(i)));
xlabel('Frequency (Hz)')
ylabel('Length of the spring (m)')
zlabel('Response Amplitude (m)')
end
Donghun Lee
2020-5-12
Ameer Hamza
2020-5-12
I am glad to be of help.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Creating, Deleting, and Querying Graphics Objects 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
