MATLAB: How to organise the code

graphplotplotting

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.

Best Answer

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
Related Question