Using ODE Event Location and using Solution at Event as Input into ode45

Question

I have been attempting to re-create a plot found in the following paper:

http://iopscience.iop.org/article/10.1088/0964-1726/4/4/006/meta

The plot I am trying to re-create is on page 267 top right hand corner, specifically, $x(t)$ vs. $t$.

The ODE used in the paper is $$m\ddot{x}(t)+k_0x(t)+k_1(x(t)-x_s)=0$$ where $x_s$ is a piece-wise constant function defined by $x_s=x(t_r)$ whenever $\dot{x}(t_r)=0$.

The plot is generated using the values $m=1,k_0=150,$ and $k_1=50.$

The code to solve this is given as:

tstart = 0;
tfinal = 2;

tspan = [tstart tfinal];
u0 = [1; 0];
refine = 4;
options = odeset('Events', @resetevent , 'OutputSel' , 1, 'Refine', refine);

tout = tstart;
uout = u0';
teout = [];
ueout = [];
ieout = [];
ue = [];
while tout(end) < tfinal

% Solve until the first terminal event.
[t, u, te, ue, ie] = ode45(@(t,u) sys(t, u, ue), [tstart tfinal] , u0, options);

% Accumulate output.  This could be passed out as output arguments.

nt = length(t);
tout = [tout; t(2:nt)];
uout = [uout; u(2:nt,:)];
teout = [teout; te];          % Events at tstart are never reported.
ueout = [ueout; ue];
ieout = [ieout; ie];


% Set the new initial conditions, with .9 attenuation (??).

u0 = [u(nt,1) ; u(nt,2)];

% A good guess of a valid first timestep is the length of the last valid
% timestep, so use it for faster computation.  'refine' is 4 by default.
options = odeset(options,'InitialStep',t(nt)-t(nt-refine),...
    'MaxStep',t(nt)-t(1));

tstart = t(nt);
end

figure(1)
set(gcf,'units','normalized','outerposition',[0 0 1 1])

subplot(2,1,1);
plot(tout,uout(:,1));
title(['Time history of $x(t)$'],'interpreter','latex','fontsize',20)
ylabel('$x$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

subplot(2,1,2);
plot(tout,uout(:,2))
title('Time history of $\dot{x}(t)$','interpreter','latex','fontsize',20)
ylabel('$\dot{x}$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

subplot(2,1,1);
stairs(teout,ueout(:,1))
title('Time history of $x_s(t)$','interpreter','latex','fontsize',20)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

% 
function dz = sys(t, z , us)
dz = zeros(2,1);
if isempty(us) == 1
    us = 0;
else
    us;
end

% if isempty(us) == 1
%     us = 0;
% elseif us(1) == 1
%     us(1) = -1;
% else
%     us;
% end

m = 1;
k0 = 150;
k1 = 50;

dz(1) = z(2);
dz(2) = - ((k0 + k1)/m)*z(1) + (k1/m)*us(1);
end


% Reset Event Function
function [lookfor, stop, direction] = resetevent(t, z)
% Locate the time when xdot passes through zero in any direction
% and stop integration.
lookfor = z(2);  % when xdot = 0
stop = 1;        % stop the integration
direction = 0;   % any direction
end

The plot produced with u0 = [1 ; 0] is

The plot produced with u0 = [1 ; -eps] is

The plot produced with u0 = [u(nt,1) ; 0] as the restart initial condition is

I have a couple questions related to this:

1. What's confuses me is that for some reason Matlab does actually find the first zero of $\dot{x}$ and inputs into the function, but the documentation on ODE Event Location says it skips the first terminal event. It's probably something I'm doing wrong, but I can't see what it is.

2. If you look at plot for $\dot{x}(t)$, at the points where it passes through zero, I have the initial conditions start with the last value it ended with and it gives a kink in the plot, which means the derivatives don't match. I tried using 0 as the initial value, which makes sense to me, but it causes the solver to diverge. I'm not sure what the 0.9 attenuation does and I tried messing with that to no avail.

If you use the commented code in the function sys it will force the value of us and the plot will look like the one in the paper, but it's too specific and not what I'm trying to do, since I want to extend this idea of resetting to a different ODE.

Thank you for any help!

Answer 1

You can work around this problem by moving forward a short duration in time, for example $10^{-9}/|\ddot{x}|$, after each event, using a call to ode45 with no event handling. This will move you forward along the curve to a point where $\dot{x}$ is definitely not zero (unless your solution is somehow quite nasty). Then you can restart there with event handling re-enabled.

In particular, since the plot in the original paper has an extremum at the very start but the "control" is not supposed to kick in until the first "interior" extremum, you will want to move forward a little bit from the initial condition before enabling event handling as well.

It's not very clear to me why this is necessary; it seems like it would be routine to find an event, move forward to outside some neighborhood of the event, and then continue searching for an event.