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.

Answer 1

[Ameer Hamza]

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

[Donghun Lee]

Sorry for the late check. I really appreciate your help Ameer. By the way, this code doesn't give us clear images of the graph for individual value of d (i.e. the graph looks like a single graph rather than the combinations of different graphs). How do I make it clear for the graph to represent 5 different graphs (in this case we have 5 values of d) in one chart? Hope you understand what I mean by. Thank you for your time.

[Ameer Hamza]

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]

Thank you so much Ameer! It was very helpful

[Ameer Hamza]

I am glad to be of help.