How do I organise my code?

clc
clear;

k_l = 26400; %Linear stiffness
m = 483; %Mass
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
dv = linspace(0,-1,40); % loop values for l
l = 0.01;
maxRes = zeros(length(dv),1);
count = 1;

for d = dv
    Om_array = linspace(0,20,10); %in rad/s-1
    A_array = linspace(0,0.06,10);
    [om_array, a_array] = meshgrid(Om_array, A_array);
    Response_amp = zeros(size(Om_array));
    T = 130;
    x0 = [0,0];
    for i=1:numel(Om_array)
        for j=1:numel(A_array)
            Om = om_array(i,j);
            A = a_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
 
    maxRes(count) = max(Response_amp(:));
    count = count+1;
end

%%%%%%%%%%%%%%figure(count)
grid on
plot(dv,maxRes, 'LineWidth', 1.5)


hold on


k_l = 26400; %Linear stiffness
m = 483; %Mass
dv = linspace(0,-1,40); % loop values for l
maxRes = zeros(length(dv),1);
count = 1;
l = 0.02;

for d = dv
    Om_array = linspace(0,20,10); %in rad/s-1
    A_array = linspace(0,0.06,10);
    [om_array, a_array] = meshgrid(Om_array, A_array);
    Response_amp = zeros(size(Om_array));
    T = 130;
    x0 = [0,0];
    for i=1:numel(Om_array)
        for j=1:numel(A_array)
            Om = om_array(i,j);
            A = a_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
 
    maxRes(count) = max(Response_amp(:));
    count = count+1;
end

grid on
plot(dv,maxRes,'LineWidth', 1.5)


hold on


k_l = 26400; %Linear stiffness
m = 483; %Mass
dv = linspace(0,-1,40); % loop values for l
maxRes = zeros(length(dv),1);
count = 1;
l = 0.03;

for d = dv
    Om_array = linspace(0,20,10); %in rad/s-1
    A_array = linspace(0,0.06,10);
    [om_array, a_array] = meshgrid(Om_array, A_array);
    Response_amp = zeros(size(Om_array));
    T = 130;
    x0 = [0,0];
    for i=1:numel(Om_array)
        for j=1:numel(A_array)
            Om = om_array(i,j);
            A = a_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
 
    maxRes(count) = max(Response_amp(:));
    count = count+1;
end


hold on


k_l = 26400; %Linear stiffness
m = 483; %Mass
dv = linspace(0,-1,40); % loop values for l
maxRes = zeros(length(dv),1);
count = 1;
l = 0.04;

for d = dv
    Om_array = linspace(0,20,10); %in rad/s-1
    A_array = linspace(0,0.06,10);
    [om_array, a_array] = meshgrid(Om_array, A_array);
    Response_amp = zeros(size(Om_array));
    T = 130;
    x0 = [0,0];
    for i=1:numel(Om_array)
        for j=1:numel(A_array)
            Om = om_array(i,j);
            A = a_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
 
    maxRes(count) = max(Response_amp(:));
    count = count+1;
end

grid on
plot(dv,maxRes,'LineWidth', 1.5)

grid on
plot(dv,maxRes,'LineWidth', 1.5)

h = legend('0.01', '0.02', '0.03', '0.04');
h.Title.String = 'l (m)';
set(h,'FontSize',18);

xlabel('Pretension (m)')
ylabel('The maximum response amplitude (m)')
set(gca,'FontSize',17)

Hi, all. This code shows a graph (pretension vs the maximum response amplitude) depending on the value of 'l'. I wish to vary l from 0 to 0.5 but it makes my code too long and therefore it takes ages to plot the graph.

Can anyone help me out on this please?

Thank you for your time.

Answer 1

[Serhii Tetora]

You can try do it with loop and function

clc; clear; close all;

k_l = 26400; %Linear stiffness
m = 483; %Mass
dv = linspace(0,-1,40); % loop values for l
l = 0:.01:.05;
for i = 1:length(l)
    maxRes(:,i) = maxResfunc(k_l,m,dv,l(i));
end

plot(dv,maxRes,'LineWidth', 1.5)
grid on
h = legend({num2str(l')},'Location','northwest','FontSize',12);
h.Title.String = 'l (m)';
xlabel('Pretension (m)','FontSize',12)
ylabel('The maximum response amplitude (m)','FontSize',12)


function maxRes = maxResfunc(k_l,m,dv,l)
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
maxRes = zeros(length(dv),1);
count = 1;
for d = dv
    Om_array = linspace(0,20,10); %in rad/s-1
    A_array = linspace(0,0.06,10);
    [om_array, a_array] = meshgrid(Om_array, A_array);
    Response_amp = zeros(size(Om_array));
    T = 130;
    x0 = [0,0];
    for i=1:numel(Om_array)
        for j=1:numel(A_array)
            Om = om_array(i,j);
            A = a_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
    maxRes(count) = max(Response_amp(:));
    count = count+1;
end
end

[donghun lee]

Thank you so much!! It is really helpful. I appreciate your effort.