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
| |
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.670393The accuracy of the computed results is what we might expect for the specified error tolerances.
After running this program, sol.xe is
0.0000 4.5168 4.5168 9.7511 9.7511 11.6704This 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.