MATLAB: How to simplify the code

3d plotsgraphplot

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.

Best Answer

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',[])
Related Question