how do I plot 3d graph?

clc
clear all

A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
%Om = 9.512/(2*pi);

k_l = 26400; %Linear stiffness
m = 483; %Mass
%d = -0.1; %Stretching condition


f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
fig = figure();
ax = axes();
hold(ax); 
% view([-53 33]);
grid on
l_array = linspace(0.1,1,100); %Excitation Amplitude
d_array = linspace(-0.1,-1,100);
T = 150;
x0 = [0,0];
xval = zeros([],1) ; 
yval = zeros([],1) ; 
for i=1:numel(l_array)
    for i=1:numel(d_array)
    l = l_array(i);
    d = d_array(i);
    k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
    
f = @(t,x) [ x(2); ...
    -(2*k_s*(x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))* ...
    (sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2) - l))/ ...
    (m*(sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2))) ];

    [t, x] = ode45(f,[100,T],x0);
    Response_amp = (max(x(:,1)) - min(x(:,1)))/2;
%    xval(i) = Om/(2*pi) ; 
    xval(i) = l;
    yval(i) = Response_amp ; 
    zval(i) = d;
    end
    % plot(Om, Response_amp, '.');
end
plot3(xval,yval,zval) ; 



xlabel('length of spring (m)')
ylabel('Response Amplitude (m)')
zlabel('pretension (m)')
set(gca,'FontSize',13)

Hi, all. I wish to plot 3d grpah (length of spring vs response amplitude vs pretension). If run this code, I still get 2d graph. How am I supposed to do?

Thank you for your time.

Answer 1

[Ameer Hamza]

It is creating 3D line, but it is not visible because you run the hold(ax) at the beginning. Add the view(3) line like this to see the 3D view

fig = figure();
ax = axes();
view(3); % <--- add this
hold(ax);

Also, note that you are using i as a variable in both loops. I think this might be a mistake. Please check according to your program logic.

[donghun lee]

Hi, thank you for your reply. If I shouldn't use i as a variable in both loops, how am I supposed to write the code? I'm quite confused about this.

[Ameer Hamza]

You can use i for out loop and j for inner loop. Use different variable names for each for loop.

[donghun lee]

Yes, if I change the variable for inner loop, it gives a differen result. However, what does this physically mean? Why do I need to make the inner variable different from the outer variable?

Thank you for your time.

[Ameer Hamza]

The meaning of nested for loop depends on what you are trying to do. In your current code, the outer for loop is selecting the elements from the d_array and d_array array. The inner for loop is not doing anything special. I suspect currently your code is giving wrong results. Can you explain the logic behind your code, maybe I can suggest the correction in your code.

[donghun lee]

I set l_array and d_array as linspace(0.1,1,100) and linspace(-0.1,-1,100), respectively.

Then, as the values of l and d change, it will give different Response_amp in the code.

I want to represent this as a 3d plot.

Thanks.

[Ameer Hamza]

Since you have two independent variables and one dependent variable, so you are trying to make a surface plot, I have refactored your code a bit. Try it

clc
clear all

A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
%Om = 9.512/(2*pi);

k_l = 26400; %Linear stiffness
m = 483; %Mass
%d = -0.1; %Stretching condition


f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
l_array = linspace(0.1,1,100); %Excitation Amplitude
d_array = linspace(-0.1,-1,100);
[L_array, D_array] = meshgrid(l_array, d_array);
Response_amp = zeros(size(L_array));
T = 150;
x0 = [0,0];
for i=1:numel(l_array)
    for j=1:numel(d_array)
        l = L_array(i,j);
        d = D_array(i,j);
        k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
        
        f = @(t,x) [ x(2); ...
            -(2*k_s*(x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))* ...
            (sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2) - l))/ ...
            (m*(sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*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

fig = figure();
ax = axes();
view(3);
hold(ax);
% view([-53 33]);
grid on
mesh(L_array,D_array,Response_amp) ;

xlabel('length of spring (m)')
ylabel('pretension (m)')
zlabel('Response Amplitude (m)')
set(gca,'FontSize',13)

[donghun lee]

WOWW!!! This is perfect! you have just saved my life! Thank you soooo much for your help sincerely!

[Ameer Hamza]

Glad to be of help.

[donghun lee]

oh yes definitely! Thanks again.

[donghun lee]

Hi, sir. Can I ask you another quick question please?

%%
clc
clear all
A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
%Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
%Om = 9.512/(2*pi);
k_l = 26400; %Linear stiffness
m = 483; %Mass
l = 0.946;
d = -0.473;
%d = -0.1; %Stretching condition
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
Om_array = linspace(0,200,100); %Excitation Amplitude
A_array = linspace(0.1,1,100);
[Om_array, A_array] = meshgrid(Om_array, A_array);
Response_amp = zeros(size(Om_array));
T = 150;
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)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))* ...
            (sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2) - l))/ ...
            (m*(sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2))) ];
        
        [t, x] = ode45(f,[100,T],x0);
        A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);
        Response_amp(i,j) = (max(x(:,1)) - min(x(:,1)))/2;
        %    xval(i) = Om/(2*pi) ;
    end
end
%% plot
fig = figure();
ax = axes();
view(3);
hold(ax);
% view([-53 33]);
grid on
mesh(Om_array,A_array,Response_amp) ;
% xlabel('length of spring (m)')
% ylabel('pretension (m)')
% zlabel('Response Amplitude (m)')
a = colorbar;
a.Label.String = 'Response Amplitude (m)';
set(gca,'FontSize',15)

This is another code that I wish to plot, however if I run this code, it keeps saying that Index in position 2 exceeds array bounds. I don't understand what it means. How would I solve this problem?

Thanks!

[Ameer Hamza]

Donghun, I my original code, I used different variable names in these three lines. First two lines have l_array and d_array on RHS. Last line have L_array and D_array. MATLAB is a case-sensitive language.

l_array = linspace(0.1,1,100); %Excitation Amplitude
d_array = linspace(-0.1,-1,100);
[L_array, D_array] = meshgrid(l_array, d_array);

In your code, you have used same name

Om_array = linspace(0,200,100); %Excitation Amplitude
A_array = linspace(0.1,1,100);
[Om_array, A_array] = meshgrid(Om_array, A_array);

Please refer to the variables name in my code to see which variables names need to start with small case and capital case inside the for loop.

[donghun lee]

Thanks Ameer. Yes, I have changed the variables deliberately because I wanted to plot a different graph, which is Om_array vs A_array vs Response_amp, not l_array vs d_array vs Response_amp. However, it just doesn't work. How am I supposed to deal with this?

Thank you for your time.

[Ameer Hamza]

Donghun, please check this code. I have modified your second code. This time I used a more robust approach

%%
clc
clear all
A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
%Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
%Om = 9.512/(2*pi);
k_l = 26400; %Linear stiffness
m = 483; %Mass
l = 0.946;
d = -0.473;
%d = -0.1; %Stretching condition
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
Om_array = linspace(0,200,100); %Excitation Amplitude
A_array = linspace(0.1,1,100);
[Om_array, A_array] = meshgrid(Om_array, A_array);
Response_amp = zeros(size(Om_array));
T = 150;
x0 = [0,0];
for i=1:size(Om_array,1)
    for j=1:size(Om_array,2)
        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)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))* ...
            (sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2) - l))/ ...
            (m*(sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2))) ];
        
        [t, x] = ode45(f,[100,T],x0);
        A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);
        Response_amp(i,j) = (max(x(:,1)) - min(x(:,1)))/2;
        %    xval(i) = Om/(2*pi) ;
    end
end
%% plot
fig = figure();
ax = axes();
view(3);
hold(ax);
% view([-53 33]);
grid on
mesh(Om_array,A_array,Response_amp) ;
% xlabel('length of spring (m)')
% ylabel('pretension (m)')
% zlabel('Response Amplitude (m)')
a = colorbar;
a.Label.String = 'Response Amplitude (m)';
set(gca,'FontSize',15)

[donghun lee]

Hi, Thank you for your reply! However, it doesn't seem like working. Can you please confirm this again if you are available?

Thank you for your time Ameer.

[Ameer Hamza]

Donghun, it gives an answer on my machine, although it takes quite a long time. Can you explain the error you get?

[donghun lee]

Hi! Sorry, it was my fault. It just takes quite long to run. If you don't mind, can I ask you a last question?

%%
clc
clear all
A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
%Om = 9.512/(2*pi);
k_l = 26400; %Linear stiffness
m = 483; %Mass
l = 0.946;
d = -0.473;
%d = -0.1; %Stretching condition
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
%Om_array = linspace(0,200,100); %Excitation Amplitude
A_array = linspace(0,1,100);
% [Om_array, A_array] = meshgrid(Om_array, A_array);
Response_amp = zeros(size(A_array));
T = 150;
x0 = [0,0];
for i=1:size(A_array)
        A = A_array(i);
        k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
        
        f = @(t,x) [ x(2); ...
            -(2*k_s*(x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))* ...
            (sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2) - l))/ ...
            (m*(sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2))) ];
        
        [t, x] = ode45(f,[100,T],x0);
        A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);
        Response_amp(i) = (max(x(:,1)) - min(x(:,1)))/2;
        %    xval(i) = Om/(2*pi) ;
    end

%% plot
fig = figure();
ax = axes();
%view(3);
hold(ax);
% view([-53 33]);
grid on
plot(A_array,Response_amp) ;
 xlabel('A (m)')
 ylabel('Response Amplitude (m)')
% zlabel('Response Amplitude (m)')
a = colorbar;
a.Label.String = 'Response Amplitude (m)';
set(gca,'FontSize',15)

In this code, I have defined that A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t). What I'm trying to do is to plot A vs Response_amp. However, this code makes Response_amp to be 0 all the way up to A = 1m. How could this happen?

Finally, I appreciate your efforts again.

Thanks, Ameer.

[Ameer Hamza]

donghun, are you ploting A vs Response_amp in this line?

plot(A_array,Response_amp);

this is not correct.

But there are few confusion about your code. Why do you need for loop, Nothing is changing inside the for loop. You are changing value of A

A = A_array(i);

But that value is not used anywhere else in the loop. Also, what is the purpose of this line

A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);

[donghun lee]

I'm trying to plot A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t) vs Response_amp.

The purpose of line below is just for the simplification.

A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);

[donghun lee]

I want to see the effects of changing A on Response_amp.

[Ameer Hamza]

But A is not used in your ode45 equation. So everytime you run the for loop, you get the same result. For loop runs 100 times but the output is same. Is your ode45 equation correct?

[donghun lee]

If I make A to be used in my ode45, the code I have generated is as below,

clc
clear all
A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
%Om = 9.512/(2*pi);
k_l = 26400; %Linear stiffness
m = 483; %Mass
l = 0.946;
d = -0.473;
%d = -0.1; %Stretching condition
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
%Om_array = linspace(0,200,100); %Excitation Amplitude
A_array = linspace(0,5,100);
% [Om_array, A_array] = meshgrid(Om_array, A_array);
Response_amp = zeros(size(A_array));
T = 150;
x0 = [0,0];
for i=1:size(A_array)
        A = A_array(i);
        k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
        
        f = @(t,x) [ x(2); ...
            -(2*k_s*(x(1)-A))* ...
            (sqrt((l-d)^2 + (x(1)-(A))^2) - l)/ ...
            (m*(sqrt((l-d)^2 + (x(1)-(A))^2))) ];
        
        [t, x] = ode45(f,[100,T],x0);
        A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);
        Response_amp(i) = (max(x(:,1)) - min(x(:,1)))/2;
        %    xval(i) = Om/(2*pi) ;
    end

%% plot
fig = figure();
ax = axes();
%view(3);
hold(ax);
% view([-53 33]);
grid on
plot(A,Response_amp) ;
 xlabel('A (m)')
 ylabel('Response Amplitude (m)')
% zlabel('Response Amplitude (m)')
a = colorbar;
a.Label.String = 'Response Amplitude (m)';
set(gca,'FontSize',15)

However, it says 'Vectors must be the same length'. I still don't get this.

[Ameer Hamza]

See this code. I changed two lines

• First line is for i=1:numel(A_array), please see the documentation of numel and size to see which one is appropriate in specific scenarios.
• plot(A_array,Response_amp), A_array is 1x100 and Response_amp is also 1x100.

clc
clear all
A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
%Om = 9.512/(2*pi);
k_l = 26400; %Linear stiffness
m = 483; %Mass
l = 0.946;
d = -0.473;
%d = -0.1; %Stretching condition
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
%Om_array = linspace(0,200,100); %Excitation Amplitude
A_array = linspace(0,5,100);
% [Om_array, A_array] = meshgrid(Om_array, A_array);
Response_amp = zeros(size(A_array));
T = 150;
x0 = [0,0];
for i=1:numel(A_array)
    A = A_array(i);
    k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
    
    f = @(t,x) [ x(2); ...
        -(2*k_s*(x(1)-A))* ...
        (sqrt((l-d)^2 + (x(1)-(A))^2) - l)/ ...
        (m*(sqrt((l-d)^2 + (x(1)-(A))^2))) ];
    
    [t, x] = ode45(f,[100,T],x0);
    A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);
    Response_amp(i) = (max(x(:,1)) - min(x(:,1)))/2;
    %    xval(i) = Om/(2*pi) ;
end

%% plot
fig = figure();
ax = axes();
%view(3);
hold(ax);
% view([-53 33]);
grid on
plot(A_array,Response_amp) ;
xlabel('A (m)')
ylabel('Response Amplitude (m)')
% zlabel('Response Amplitude (m)')
a = colorbar;
a.Label.String = 'Response Amplitude (m)';
set(gca,'FontSize',15)

[donghun lee]

Ahh.. Thank you sooo sooo much for your effort. This is perfect. I don't know what to say more. I appreciate you sincerely Ameer.

Best wishes.

[Ameer Hamza]

Glad to be of help.