Problem with the "Plot" function

Hi, I need some help from those of you who are good at Matlab sinc I am not a profession.. I appreciate the time and effort.

I need to plot my time in x axis and my roll angle(x1) as y axis

----This is my parameters------- function dx = wheelchaircorrectversion (t,x) global m g d F h Ix w m=10; g=-9.81; F=5; d=10; h=20; Ix=30; w=20;

if t>0 && t>=0.01 F=5 elseif t>0.01 F=0 end

z=(sqrt((.5*w).^2+(.5*h).^2)).*(sin(atan((h/2)/w/2))); d=(sqrt((.5*w).^2+(.5*h).^2)).*(cos(atan((h/2)/w/2)));

dx=zeros(2,1); dx(1)=x(2); dx(2)=(-(m*g*d)+(F+z))/Ix; end

---- my main function---

clear all clc global m g d F h Ix w m=10; g=-9.81; F=5; d=10; h=20; Ix=30; w=20;

x1_ini = [0 0]'; t_dim = 1; i=1; dt = 0.01; duration = 0.05;

for t = 0.0 : dt : duration t0 = t; i=i+1; tf = (t+dt)*(1+eps); tspan = [t0 tf];

    [t1,x1] = ode45('wheelchaircorrectversion',tspan,x1_ini);
  n1 = length(x1); %length of a vector
  tmp1 = x1(n1,:); %the last value
    tmp2 = t1(t_dim)

x1_ini = tmp1; % update initial conditions

    traj_p0(i,:) = x1_ini;
    traj_t(i)= tmp2
end

subplot(211) %Creates axes in tiled positions plot(traj_t,x1(:,1)) xlabel('Time') ylabel('Roll Angle')

subplot(212) plot(traj_t,x1(:,2)) xlabel('Time') ylabel('Roll Velocity')

----

Whenever I try to plot these functions, it says that vectors needs to be in same length for plot(traj_t, x1(;,1)) and for plot(traj_t,x1(:,2))

I have tried everything to make the number of vectors equal to each other, but it didnt work..

could someone help?

[Jan]

Please format the code properly. Thanks.

Answer 1

[Image Analyst]

More accurately, x1(:,1) needs to be the same number of elements at traj_t. x1(:,2) needs to be that length also. Do this before you call plot

size(traj_t)
size(x1)

Tell us what it says.

[Tae Yeong Kim]

Error using plot Vectors must be the same lengths.

Error in wheelchaircvparameters (line 39) plot(traj_t,x1(:,1))

I am still getting the same error.

[Image Analyst]

Repeating:

 Do this before you call plot

size(traj_t)
size(x1)

Tell us what it says. There will be some numbers in the command window that are the sizes. You could also try this:

whos traj_t
whos x1

Do that too and tell us what it says.

[Tae Yeong Kim]

tmp2 =

    0.0500

traj_t =

         0         0    0.0100    0.0200    0.0300    0.0400    0.0500

ans =

     1     7

ans =

    41     2

Name        Size            Bytes  Class     Attributes

traj_t      1x7                56  double              

Name       Size            Bytes  Class     Attributes

x1        41x2               656  double

This is what it shows me after ridiculously large amount of F=5

[Image Analyst]

OK, and you're plotting 7 values of traj_t against 42 values of x1. Does that make sense to you, because it doesn't to me?

[Tae Yeong Kim]

it doesn't make any sense.. I found n, which is the number of the vectors, and it appeared to be 42.

I tried to make the number of values of traj_t as 42, but I did not know how to make it..

I appreciate your time and effor very much

[Image Analyst]

Set a breakpoint and type this into the command window:

t = 0.0 : dt : duration

See how many elements it is. That is the loop where you're setting traj_t so that is how many elements traj_t will have. If you don't like it, you have to adjust dt or duration.

[Tae Yeong Kim]

Wow, thanks for your help. I finalized it, and it is working. However, I'm not sure if it's physically making sense. Do you have any thoughts about my code?

[Image Analyst]

Nope - I didn't even run it.

Answer 2

[Youssef Khmou]

Try this :

1.Delete the traj_t variable that exists inside the loop.

2. create a vector time, outside the loop, of length 41 such as : time=linspace(0,duration,41).

3. replace,in the plot section, traj_t wit time'.

The code now works well, but are the results reasonable In Physic viewpoint ?

As you posted non organized code , i post the new version of the code here for further discussions by other users :

--------------------------------------------------------------------------

function dx = wheelchaircorrectversion (t,x)

global m g d F h Ix w

m=10; g=-9.81; F=5; d=10; h=20; Ix=30; w=20;

if t>0 && t>=0.01

    F=5;

elseif t>0.01

    F=0;

end

z=(sqrt((.5*w).^2+(.5*h).^2)).*(sin(atan((h/2)/w/2)));

d=(sqrt((.5*w).^2+(.5*h).^2)).*(cos(atan((h/2)/w/2)));

dx=zeros(2,1); dx(1)=x(2); dx(2)=(-(m*g*d)+(F+z))/Ix;

end -------------------------------------------------------------------

And the main program :

-------------------------------------------------------------------

clear all ,clc

global m g d F h Ix w

m=10;g=-9.81; F=5; d=10; h=20; Ix=30; w=20;

x1_ini = [0 0]'; t_dim = 1; i=1; dt = 0.01; duration = 0.05;

for t = 0 : dt : duration

    t0 = t;i=i+1; tf = (t+dt)*(1+eps); tspan = [t0 tf];

    [t1,x1] = ode45('wheelchaircorrectversion',tspan,x1_ini);

n1 = length(x1); %length of a vector

tmp1 = x1(n1,:); %the last value

    tmp2 = t1(t_dim);

    x1_ini = tmp1; % update initial conditions

    traj_p0(i,:) = x1_ini;

    %traj_t(i)= tmp2;

end

time=linspace(0,duration,41);time=time';

subplot(211), plot(time,x1(:,1)), xlabel('Time'), ylabel('Roll Angle')

subplot(212), plot(time,x1(:,2)), xlabel('Time'), ylabel('Roll Velocity')

--------------------------------------------------------------------------

[Tae Yeong Kim]

Thank you for your response. I have one more question about my code: for the last part,subplot(211), plot(time,x1(:,1)), xlabel('Time'), ylabel('Roll Angle'), I am trying to plot the function before it went through ODE45 so that I get un-integrated (original) function to get roll velocity on subplot 212.

I am not sure how I have to write out the function..

Thanks a lot!

Answer 3

[Youssef Khmou]

Tae Yeong Kim,

so its about ordinary differential equation ! THE TIME IS ALREADY THERE !!

[t1,x1]=ode('wheelchaircorrectversion',tspan,x1_ini);

Position and velocity are in the 41x2 matrix x1 and the time axis is t1(41,1)

you can replace the vector "time" that we created with t1, ITS THE SAME .

About your last question i do not understand what you mean : you have a position X and Velocity V , you said you want to integrate the x1(:,1) to get the Velocity, i think its wrong :

YOU HAVE POSITION X , VELOCITY V and ACCELERATION A

DX/Dt= V(t) , DV/Dt= A, INTEG(V(t))=X(t), INTEG(A)= V(t)

so lets see :

---------------------------------------------------------------------------

plot(t1,x1(:,2)) % This is the velocity

hold on

Velocity2=diff(x(:,1))./diff(t1); % Diff function truncates to n-1

Velocity2(end+1)=Velocity2(end); % add the last element

plot(t1,Veclocity2,'r'),

-----------------------------------------------------------------------

Numerically its the same .

Advice :

1.do not declare the variables again the testing script .

2. try to comment every variable to explain to the reader , for example :

------------------------------------------------------------------------

m=10; % mass in Kg

g=-9.81; % Gravity acc in m/s²

F=5; % Force in Newton

d=10; % Distance between ..and .. in meteres....

I hope that helps

KHMOU Youssef,

[Tae Yeong Kim]

Thank you so much. It answered everything I needed.

If you could check out my equations and see if they make sense, it would be awesome.

I am not sure if variable z and d are correctly derived.. Thanks

Answer 4

[Youssef Khmou]

Hi

i think there is a problem, based on what you wrote :

Z/d=h/w. and based on numerical values u gave h=w=20. in that case z=d which is wrong ,

what is the point in the center of the rectangle? what is its dynamic ?

[Tae Yeong Kim]

That is true.. I didn't see it

If h and w were different, would the equation just be zw = dh?

The point in the center of the rectangle is the center of mass, and I have not defined it yet; if I were to define it, I would use (0,0)

I drew the free diagram to show the roll angle and roll velocity for this object to start an analysis of wheelchair dynamics. I know it is a lot basic, but I will add springs and dampers in the future.

If you have any ideas for it, it would be awesome!

Thanks

[Tae Yeong Kim]

Oh, and I forgot to tell you that z and d changes as it the object is moving from the effect of the force.

So they are variables, but I do not know the right equations for them if w and h are not equal to each other.

[Jan]

Please, Tae Yeong Kim, take the time to format your code. Currently it looks very messy.

[Tae Yeong Kim]

clear all 
clc 
global m g F h Ix w 
m=10; 
g=9.81; 
F=5; 
h=20; 
Ix=30; 
w=20;

x1_ini = [0 0]';
i=1; 
dt = 0.01; 
duration = 0.05;

for t = 0.0 : dt : duration 
    t0 = t;
    i=i+1;
    tf = (t+dt)*(1+eps);
    tspan = [t0 tf];

      [t1,x1] = ode45('wheelchaircorrectversion',tspan,x1_ini);
    n1 = length(x1); %length of a vector
    tmp1 = x1(n1,:); %the last value

      n2 = length(t1);
      tmp2 = t1(n2);

      x1_ini = tmp1; % update initial conditions

      traj_p0(i,:) = x1_ini;
      traj_t(i)= tmp2;
  end

timetraj = traj_t';
timetraj = linspace(0,duration,41)% Generates "duration" points between X1 and X2

subplot(211) %Creates axes in tiled positions 
plot(timetraj,x1(:,1)) 
xlabel('Time') 
ylabel('Roll Angle')

subplot(212) 
plot(timetraj,x1(:,2)) 
xlabel('Time') 
ylabel('Roll Velocity')
----------------------------------------------------------------------

Does this look any better?