MATLAB Answers

Tae Yeong Kim

How to use ODE 45 to simulate a model with two state equation?

Asked by Tae Yeong Kim
on 13 Jan 2013

Hi, I have started to use MATLAB from 2 months ago, and I am having a trouble of simulating a box model with parameters through ODE 45.

I cannot get the right plots of graph, even though I tried so many times to debug my code.

My following is the code that I set up, however, I do not know what is wrong from my codes

------- Here is the variables--

function dx = wheelchaircorrectversion (t,x)
global m g d F h Ix w
if t==0.01
else t>=0.02 

---- and this is the main function

clear all
global m g d F h Ix w
x1_ini = [0 0]';
dt = 0.01;
duration = 0;
for t = 0.0 : dt : duration
	i = i+1;
	t0 = t;
	tf = (t+dt)*(1+eps);
    tspan = [t0 tf];
      [t1,x1] = ode45('wheelchaircorrectversion',tspan,x1_ini);
  	n1 = length(x1);
  	tmp1 = x1(n1,:);
	x1_ini = tmp1; % update initial conditions 
      traj_p0(i+1,:) = x1_ini;
xlabel( 'disp. in x direction(m)' )
ylabel( 'disp. in y direction(m)' )


As you can see above, the codes are representing a simple box model, pushed with a force of 5 Newtons to simulate the roll angle of the object: x1 being phi, and x2 being phi dot.

The parameters are randomly put, and the time intervals are set so that I could get the control over time to put force for only a 0.0001 second, not continuously.

I am a beginner in MATLAB so I have only little knowledge of coding.. Could anyone see the problem I am having?

Thank you for your time.

[EDITED, Jan, code formatted]



No products are associated with this question.

1 Answer

Answer by Jan Simon
on 13 Jan 2013
Edited by Jan Simon
on 13 Jan 2013

ODE45, as the other ODE integrators of Matlab, requires a continuously differentiable function. Otherwise the stepsize control can explode. You find several exhaustive explanation in this forum for this topic, but sometimes I'm afraid I'm the only one who cares about this detail.

This looks strange:

if t==0.01
else t>=0.02 

Do you mean "elseif"?

Using global variables to define parameters is suboptimal. See benefit of using anonymous functions .


Yes, I thought the one time loop is necessary, because the force needed to be applied only once, not as a continuous force.

So basically, I am trying to simulate a 10 kg wheelchair, symbolized as a box, with 5 N force.

The following is the free body diagram and the variable that I have set for the diagram. Also , I set

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)));

as variables that changes as the object rotates.

Did this help you to understand it a little better?

Jan Simon
on 14 Jan 2013

No, I cannot understand it better. I cannot recognize, where try to apply a force for 0.0001 seconds. Creating a loop, which runs for 1 iteration also works, but is meaningless. I cannot see, how such a loop should support, that the force is not applied continuously. I still do not see any reaction to my argument, that the function to be integrated must be smooth or that anonymous functions have benefits.

You set i=1, then i=i1 inside the loop and write to traj_p0(i+1,:) finally. Then the first 2 rows of traj_p0 are undefined. Is this wanted?

What is "z", "d" and "w" in your formula?

Why is the "wheelchair" a square? Do you consider the rotational energy of the wheel? If not, simulation the center of mass as a sliding point would be much easier and you can even get the results in clodes form without a numerical integration.

Because I am a beginner in MATLAB, I decided to figure the wheelchair as a box that represent the main body of the wheelchair, and my finalized goal is to optimize the roll(from roll, pitch, yaw motion) stability of it. I do not know how to get the results without the integration..

Also, I know creating the loop may be meaningless, but I am doing it to get the control over the time step to 0.001 second. I have got some help from my professor for the traj_p0(i+1;) part and I think that is one of the errors that I need to fix: define the traj_p0.

w is the original width and it is constant. d represents the distance between the tipping point and the center of mass that changes overtime when force is implied, and z represents the height from ceter of mass to the contacting face.

I am sorry that I can't explain the codes thoroughly.. I am in a process or learning matlab. It would be awesome if you have any ideas to the project. Thanks again

Discover MakerZone

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

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

MATLAB Academy

New to MATLAB?

Learn MATLAB today!