MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today**New to MATLAB?**

Asked by Tae Yeong Kim
on 27 Jan 2013

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?

*No products are associated with this question.*

Answer by Image Analyst
on 27 Jan 2013

Edited by Image Analyst
on 27 Jan 2013

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.

Show 5 older comments

Image Analyst
on 27 Jan 2013

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
on 27 Jan 2013

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?

Answer by Youssef Khmou
on 27 Jan 2013

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
on 27 Jan 2013

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 by Youssef Khmou
on 28 Jan 2013

Edited by Youssef Khmou
on 28 Jan 2013

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
on 28 Jan 2013

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 by Youssef Khmou
on 28 Jan 2013

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 ?

Show 1 older comment

Tae Yeong Kim
on 29 Jan 2013

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 Simon
on 29 Jan 2013

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

Tae Yeong Kim
on 2 Feb 2013

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?

## 1 Comment

## Jan Simon (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/59991#comment_125491

Please format the code properly. Thanks.