Example 7

For some problems the changes in the equations occur at times that are not known in advance. The event location capability is used to determine when there is a change and the integration is terminated. It is then restarted with the new definition of the equation. The role of the history and the possibility of a jump discontinuity in the solution itself complicate this, but dde23 was designed to make it as painless as possible. This example and the next show how to proceed.

Marriott and DeLisle [10] solve a DDE that involves a step function of the history term. With D = y(t-12) - xb , the equation is
y(t) = ( - y(t) + p( a + e  sign(D) - u  sin2 (D) ) )/t .
It is solved on [0,120] with history y(t) = 0.6 for t 0 and parameter values xb = -0.427, a = 0.16, e = 0.02, u = 0.5,t = 1 .

The term sign(D) is an idealization of a sharp change that we do not expect to resolve in detail. dde23 is sufficiently robust that it can solve this problem in a straightforward way. However, ignoring discontinuities can get you into trouble even when solving an ODE [19], much less a DDE. You can be much more confident of your numerical results if you arrange that the solver is always integrating an equation with smooth coefficients. This can be done here by letting the parameter state be the value the solver is to use for sign(D). With the given history, it is initialized to +1. We integrate until the event function y(t-12) - xb vanishes. The 'Events' option is used to locate this event and terminate the integration then. The sign of state is changed and the integration restarted. This continues until the end of the interval is reached. It is straightforward to code the event location as

 function [value,isterminal,direction] = exam7e(t,y,Z,state)
 xb = -0.427;
 value = Z - xb;
 isterminal = 1;
 direction = 0;
and with state defined in exam7.m as described, the evaluation of the DDE is coded as
 function v = exam7f(t,y,Z,state)
 xb = -0.427;
 a = 0.16;
 epsilon = 0.02;
 u = 0.5;
 tau = 1;
 Delta = Z - xb;
 v = (-y + pi * (a + epsilon*state - u*sin(Delta)^2)) / tau;

The event location capability deals with the issue of finding when the equations change, but there is a matter special to DDEs on restarting, the history. On a restart, dde23 accepts the previously computed solution structure as history. dde23 updates the information in the solution structure each time it is called, so the solution returned is always valid from the initial to the last point reached in the integration, namely sol.x(end). It is convenient to compare this point to the end of the interval of integration and perform the restarts in a while loop. In this loop the solution structure of one integration is used as the history for the next until the run is completed.

With exam7e.m and exam7f.m, the complete program exam7.m to compute and plot the solution is

 state = +1;
 opts = ddeset('Events','exam7e');
 sol = dde23('exam7f',12,0.6,[0, 120],opts,state);
 while sol.x(end) < 120
   fprintf('Restart at %5.1f.\n',sol.x(end));
   state = - state;
   sol = dde23('exam7f',12,sol,[sol.x(end), 120],opts,state);
 end

 plot(sol.x,sol.y);
 title('Figure 7. Marriott-DeLisle Controller Problem')
 xlabel('Restart each time y(t - 12) = -0.427.');
The program prints out a message about restarts and where they occur. The resulting output is
 Restart at  14.0.
 Restart at  24.9.
 Restart at  68.1.
 Restart at  75.7.
 Restart at 118.9.

images/7.gif

References

[10]
C. Marriott and C. DeLisle, Effects of discontinuities in the behavior of a delay differential equation, Physica D, 36 (1989) 198-206.

[19]
L.F. Shampine and S. Thompson, Event location for ordinary differential equations, Comp. & Maths. with Appls., to appear.