Example 4

It is often necessary to find when a solution satisfies a certain relation, e.g., when a component has a specific value. An event is said to occur when a function of the solution, g(t,y(t),y(t -t1),,y(t - tk)), vanishes. Some problems involve many of these ``event functions''. This example shows how to use the powerful event location capability of dde23.

Figure 15.6 of [6] displays the solution of an infectious disease model. The equations
y1(x)
=
- y1(x)y2(x-1) + y2(x-10)
y2(x)
=
y1(x)y2(x-1) - y2(x)
y3(x)
=
y2(x) - y2(x-10)
are solved on [0,40] with history y1(x) = 5,y2(x) = 0.1,y3(x) = 1 for x 0. To illustrate event location, we compute the local maxima of all three solution components.

We compute the maxima by finding where the first derivatives vanish. The three event functions come from the DDEs: y1(x) = -y1(x)y2(x-1)+y2(x-10), and so forth. All event functions are evaluated in a single MATLAB function that returns the values as a column vector. The name of this function is passed to the solver as the value of the 'Events' option. For this example we evaluate the three functions in exam4e by a call to exam4f. Because event location is used for a variety of purposes, we have to tell dde23 more about what we want to do. Sometimes we just want to know that an event has occurred and other times we want to terminate the integration then. We tell the solver about this by returning a vector isterminal from exam4e. To terminate the integration when event function k vanishes, we set component k of isterminal to 1 (true), and otherwise to 0 (false). For this example none of the events is terminal. There is an annoying matter of some importance: Sometimes we want to start an integration with an event function that vanishes at the initial point. Imagine, for example, that we fire a model rocket into the air and we want to know when it hits the ground. It is natural to use the height of the rocket as a terminal event function, but it vanishes at the initial time as well as the final time. dde23 treats an event at the initial point in a special way. The solver locates such an event and reports it, but does not treat it as terminal, no matter how isterminal is set. The example shows that how an event function vanishes may be important: To distinguish maxima from minima, we want the solver to report that a derivative vanished only when it changes from positive to negative values. This is done using direction. If we are interested only in events for which event function k is increasing through 0, we set component k of direction to +1. Correspondingly, we set it to -1 if we are interested only in those events for which the event function is decreasing, and 0 if we are interested in all events. Once we understand what information must be provided, it is easy to code the event functions of this example as

 function [value,isterminal,direction] = exam4e(x,y,Z)
 value = exam4f(x,y,Z);
 isterminal = zeros(3,1);
 direction = -ones(3,1);

Now that we have discussed how to tell the solver what we want it to do, we have to discuss how it reports what happened. The locations of events are returned as the field sol.xe and the values of the solution at these points are returned as the field sol.ye. If there are no events, sol.xe = []. The field sol.ie reports which event occurred. A value of k indicates that event function k vanished at the corresponding entry of sol.xe.

It is straightforward to code the equations as

 function v = exam4f(x,y,Z)
 ylag1 = Z(:,1);
 ylag2 = Z(:,2);
 v = zeros(3,1);

 v(1) = -y(1)*ylag1(2) + ylag2(2);
 v(2) =  y(1)*ylag1(2) - y(2);
 v(3) =  y(2) - ylag2(2);
With exam4e.m and exam4f, it is also straightforward to code the solution of the problem as the first two lines of the complete solution exam4.m that follows:
 options = ddeset('Events','exam4e');
 sol = dde23('exam4f',[1, 10],[5; 0.1; 1],[0, 40],options);

 xe = sol.xe;
 ye = sol.ye;
 ie = sol.ie;

 n1 = find(ie == 1);
 x1 = xe(n1);
 y1 = ye(1,n1);
 n2 = find(ie == 2);
 x2 = xe(n2);
 y2 = ye(2,n2);
 n3 = find(ie == 3);
 x3 = xe(n3);
 y3 = ye(3,n3);

 plot(sol.x,sol.y,'k',x1,y1,'rs',x2,y2,'rs',x3,y3,'rs')
 title('Figure 4. Infectious disease model from Hairer et al.')
 xlabel('Maxima are indicated by red squares.')
The only complication in this program is separating the various kinds of events. It is not necessary, but perhaps clearer, to introduce local variables for the fields that return the results of the event location. The command n1 = find(ie == 1) finds the indices corresponding to the first event function. These indices allow us to extract the information that y1(x) has its maxima at xe(n1) and its values there are ye(1,n1). The second and third event functions are handled in the same way and then all the results are plotted.

images/4.gif

Reference

[6]
E. Hairer, S.P. Nø rsett, and G. Wanner, Solving Ordinary Differential Equations I, Springer-Verlag, Berlin, 1987.