11.3: ODE Events

Normally when you call ode45 you specify a start time and an end time. But sometimes you don’t know ahead of time when the simulation should end. To solve this problem we can define an event, something of interest that happens during a simulation, like the penny reaching the ground.

Here are the steps:

1. First we define an event function that allows ode45 to figure out when an event occurs. Here’s an event function for the penny example:

   function [value, isterminal, direction] = event_func(t,X)
       value = X(1);
       isterminal = 1;
       direction = -1;
   end

   The event function takes the same input variables as the rate function and returns three output variables: value determines when an event can occur, direction determines whether it does, and isterminal determines what happens. More specifically, an event can occur when value passes through 0. If direction is positive, the event only occurs if value is increasing. If direction is negative, the event only occurs if value is decreasing. If direction is 0, the event always occurs. If isterminal is 1, the event causes the simulation to end; if it is 0, the simulation continues.

   This event function uses the altitude of the penny as value so an event can occur when altitude passes through 0. Because direction is negative, an event occurs only when altitude is decreasing, and because isterminal is 1, the simulation ends if an event occurs.

2. Next, we use odeset to create an object called options:

   options = odeset('Events', @event_func);

   The name of the option is Events and the value is the handle of the event function.

3. Finally, we pass options as a fourth argument to ode45:

   [T, M] = ode45(@rate_func, tspan, X, options);

   When ode45 runs, it invokes event_func after each time step. If the event function indicates that a terminal event occurred, ode45 stops the simulation.

Let’s look at the results from the penny example:

>> T(end)
8.8179

>> M(end, :)
0.0000  -86.4153

The last value of T is about 8.8, which is the number of seconds the penny takes to reach the sidewalk.

The last row of M indicates that the final altitude is 0, which is what we wanted, and the final velocity is about -86 m/s.