Example 8

This example is much like Example 7 except that the solution itself is discontinuous. We restart at discontinuities, so the jump in the solution occurs at the initial point of an integration and can be handled as in Example 6.

A two-wheeled suitcase may begin to rock from side to side as it is pulled. When this happens, the person pulling it attempts to return it to the vertical by applying a restoring moment to the handle. There is a delay in this response that can affect significantly the stability of the motion. This is modeled by Suherman et alia [21] with the DDE
q(t) + sign(q(t))gcos(q(t)) - sin( q(t)) + bq(t - t) = A sin(Wt + h)
The equation is solved on [0, 12] as a pair of first order equations with y1(t) = q(t), y2(t) = q(t). Figure 3 of [21] shows a plot of y1(t) against t and a plot of y2(t) against y1(t) when g = 2.48, b = 1,t = 0.1, A = 0.75, W = 1.37, h = arcsin(g/A) and the initial history is the constant vector zero.

A wheel hits the ground (the suitcase is vertical) when y1(t) = 0 and the suitcase has fallen over when |y1(t)| = p/2. The events are terminal and all are to be reported. The event function can be coded as

 function [value,isterminal,direction] = exam8e(t,y,Z,state)
 value = [y(1); abs(y(1))-pi/2];
 isterminal = [1; 1];
 direction = [0; 0];

As in Example 7, the parameter state seen in the event function is used in exam8.m to evaluate properly the discontinuous coefficient sign(y1(t)) in the DDE. We initialize it to +1 and change its sign when dde23 returns because y1(t) vanished. However, there are two event functions, so we must check the last entry in sol.ie to see if we should change the sign of state. With this, the DDEs can be coded as

 function v = exam8f(t,y,Z,state)
 ylag = Z(:,1);
 v = zeros(2,1);

 gamma = 0.248;
 beta  = 1;
 A = 0.75;
 omega = 1.37;
 eta = asin(gamma/A);

 v(1) = y(2);
 v(2) = sin(y(1)) - state*gamma*cos(y(1)) - beta*ylag(1) ...
       + A*sin(omega*t + eta);

When a wheel hits the ground, the integration is to be restarted with y1(t) = 0 and y2(t) multiplied by the coefficient of restitution 0.913. The 'InitialY' option is used for this purpose. The solution at all the mesh points is available as the field sol.y and in particular, the solution at the time of the event is the last column of this array, sol.y(:,end). If the suitcase falls over, the run is terminated, so again we must check which event occurred. With exam8e.m and exam8f.m, the complete program exam8.m to solve the problem and plot the solution in the phase plane is

 state = +1;
 opts = ddeset('RelTol',1e-5,'AbsTol',1e-5,'Events','exam8e');
 sol = dde23('exam8f',0.1,[0; 0],[0 12],opts,state);

 % Reference values for event locations:
 ref = [4.516757065, 9.751053145, 11.670393497];
 fprintf('\nKind of Event:                 dde23   reference\n');
 event = 0;
 while sol.x(end) < 12
   event = event + 1;
   if sol.ie(end) == 1
      fprintf('A wheel hit the ground.   %10.4f  %10.6f\n',...
              sol.x(end),ref(event));
      state = - state;
      opts = ddeset(opts,'InitialY',[ 0; 0.913*sol.y(2,end)]);
      sol = dde23('exam8f',0.1,sol,[sol.x(end) 12],opts,state);
   else
      fprintf('The suitcase fell over.   %10.4f  %10.6f\n',...
              sol.x(end),ref(event));
      break;
   end
 end
 plot(sol.y(1,:),sol.y(2,:))
 title('Figure 8.  Suitcase problem.')
 xlabel('\theta(t)')
 ylabel('\theta''(t)')
Note that ddeset can be used to change the value of an option or add an option, just as with odeset. The program reproduces the phase plane plot of Figure 3 in [21]. It also reports what kind of event occurred and the location of the event. Reference values were computed with the FORTRAN 77 code DKLAG5 [13] used in [21] and verified with its successor DKLAG6 [4]. Having written the three solvers, we can fairly say that it is very much easier to solve this problem in MATLAB with dde23. The program results in the output
 Kind of Event:                 dde23   reference
 A wheel hit the ground.       4.5168    4.516757
 A wheel hit the ground.       9.7511    9.751053
 The suitcase fell over.      11.6704   11.670393
The accuracy of the computed results is what we might expect for the specified error tolerances.

images/8.gif
After running this program, sol.xe is

 0.0000    4.5168    4.5168    9.7511    9.7511   11.6704
This does not seem to agree with the event locations reported by the program. For instance, why is there an event at 0? That is because one of the event functions is y1 and this component of the solution has initial value 0. As explained in Example 4, dde23 locates this event, but does not terminate the integration because the terminal event occurs at the initial point. The first integration terminates at the first point after the initial point where y1(t*) = 0, namely t* = 4.5168. The second appearance of 4.5168 in sol.xe is the same event at the initial point of the second integration. The same thing happens at 9.7511 and finally the event at 11.6704 tells us that the suitcase fell over and we are finished.

References

[4]
S.P. Corwin, D. Sarafyan, and S. Thompson, DKLAG6: A code based on continuously imbedded sixth order Runge-Kutta methods for the solution of state dependent functional differential equations, Appl. Num. Math., 24 (1997) 319-333.

[13]
K.W. Neves and S. Thompson, Software for the numerical solution of systems of functional differential equations with state dependent delays, Appl. Num. Math., 9 (1992), 385-401.

[21]
S. Suherman, R.H. Plaut, L.T. Watson, and S. Thompson, Effect of human response time on rocking instability of a two-wheeled suitcase, J. of Sound and Vibration, 207 (1997) 617-625.