How to add two constraints in the event function of Matlab ODE's

Hello to everyone, I am trying to solve a system of diff. equations in matlab where the system must be stopped once the solution reaches two extreme values. I must do it for three bodies of the system. In order to make the things more simple for you, my problem may be compared to the ballode example implemented in matlab, where once the ball reaches the floor, the integration is stopped and restarted with new initial onditions. I am trying to add another constraint to the ball,for example a roof at height 10, in order to impede the lift further, under theese conditions, the ball should bounce back to the floor and so on until the motion has stopped. Could anyone help me to write the code ballode with this new constraint ? I suppose that new initial conditions must be added when the ball bounces from the roof and until now I havent had success.

Any help would be appreciated. Thank you very much Erjon

Answer 1

[erjon selmani]

Hello Steven, I will explain the problem in more detail and with two figures. Figure 1:

I want to see the motion (position) of an object inside a box, when some external forces are applied. If the grey square is sit down due to the resultant of these forces, I must apply the following conditions: speed = 0, acceleration = 0. Hence, the two time integration of the acceleration should give me the position (h) which is expected to be zero, until the sum of forces change and the square lifts. In Figure 2 is given the other possible condition.

When the sum of forces in upward direction overcome that in downward direction, the square is lifted up to the ceiling of the box and the motion is no further allowed due to physical reason. In this case I need to apply again the initial conditions: speed = 0, acceleration = 0, as long as the sum of forces remain the same.

As you can see, there is a certain similarity with the ballode example, but in my case I need to account a limitation of the motion in both directions while in the ballode the limit was only the floor.

Thank you

[Steven Lord]

When the height of the bottom edge of the object crosses 0 in a decreasing direction (hits the bottom of the box) event #1 triggers. Stop the solver, set the initial speed and acceleration, and restart.

When the distance between the top of the object and the top of the box crosses 0 in a decreasing direction (the object hits the top of the box) event #2 triggers. Stop the solver, set the initial speed and acceleration, and restart.

Both events are terminal. Both need to compare the position of the object with a different value. You can frame both as zero crossings from above.

There's another way to frame event #2, watching for the height of the object minus the height of the box to cross 0 from below, but then you've got two different events that watch for zero crossings in different directions. This doesn't cause a problem for the solver, but you should choose which one you think will seem more intuitive for the person reading and trying to understand the code (which could be you, six months from now.)

Answer 2

[Steven Lord]

You can specify multiple event conditions in the events function. The ballode function is the "Simple Event Location" example on this documentation page. The "Advanced Event Location" example orbitode on that same page sounds like it's closer to your actual application, and it shows how to specify multiple event conditions.

[erjon selmani]

Ok, thank you very much for the answer, i will study the orbitode matlab file.

Answer 3

[erjon selmani]

Dear @Steven Lord, I saw the orbitode matlab file but it is not what I was intending in my question. In orbitode are set two events, and the integration is stopped in correspondence to one of them, without restarting anymore. In my case I need to find when two possible events occur, to stop the integration and to restart with new initial conditions, as in the ballode file. The only difference is that I must apply two different initial conditions, say one per event, and this is the difficulty that had arisen ! Any suggestions ?

Thank you very much for the support.

Erjon

[Steven Lord]

So put the two examples together. Make a copy of ballode and modify the copy's events function to use multiple conditions like the orbitode example does.

By "I must apply two different initial conditions, say one per event" do you mean you need to set different initial conditions for the next restart based on which of your events caused the previous ODE solver call to terminate? So as an example if you were doing predator - prey modeling and the predator population goes to 0 maybe humans introduce a different predator who interacts with the prey differently, while if prey goes to 0 the predators start migrating away from the region being modeled to find food elsewhere?

If that's the case, see the "Event Information" section in the documentation page to which I linked in my answer. The fifth output of the ODE solver or the ie field of the struct returned if you call the ODE solver with one output will contain information about which event condition was detected. Use that to determine what your new initial conditions for your restart should be.

If that's not what you meant, I don't understand what you're trying to do. Can you explain in a little more detail the problem that you're trying to solve?

Answer 4

[erjon selmani]

I tried to merge your answer with the ballode file, please have a look and give your opinion.

    % code
 tstart ;
tfinal ;
h0 ;
options = odeset('Events',@events);
for i = 1:k
   % Solve until the first terminal event.
   [t,h,te,ye,ie] = ode23(@fun,[tstart tfinal],h0,options);
    if ~ishold
      hold on
    end
      h0(1) = 0;
      h0(2) = 0;
  tstart = t(nt);
   end    

 function [value,isterminal,direction] = events(t,h)
% Locate the time when height passes through zero in a decreasing direction
% and stop integration.
value = [h(1); h(1) - hmax];    % detect when the height = 0
isterminal = [1;1];   % stop the integration
direction = [1;-1];   % define the direction
% code

[erjon selmani]

The direction for the zero crossing is to be decided according to the direction of the reference system.