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

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
m=10;
g=-9.81;
F=5;
d=10;
h=20;
Ix=30;
w=20;

if t==0.01
        F=5;
else t>=0.02 
        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 this is the 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;

i=1;
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;
  end

subplot(221)
plot(traj_p0(:,1),traj_p0(:,2))
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]

Answer 1

[Jan]

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
      F=5;
else t>=0.02 
      F=0;
end

Do you mean "elseif"?

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

[Tae Yeong Kim]

Yes, I did the help function on elseif, and I think that's what I was looking for.

After changing that, the simulation is still the same.

Do you have any suggestions to revise my main function?

Again, thank you very much.. MATLAB sure is a hard apparatus to use.

[Jan]

Sorry, I cannot give further suggestions, because I do not understand, what you want to achieve. On one hand I cannot follow your test explanations, on the other hand I do not understand the code also, e.g.:

duration = 0;
for t = 0.0 : dt : duration
  ...

This runs the loop once only, but is this wanted?

As I have said already, I assume that this is not conform to the specifications of the ODE algorithm:

 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.

[Tae Yeong Kim]

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]

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.

[Tae Yeong Kim]

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