Setting ODE multiple integration start and end times whilst wanting specific intermediate points

Hi,

Until recently my model was integrated as a single time span and used a long vector of times to return values at the desired intermediate points:

t0 = 0; % s. Input signal start time
t1 = 5; % s. Input signal end time
dt = 0.005; % time resolution
tvec = t0:dt:t1; % creates a horizontal vector between t0 and t1 that increments by dt

[T, x] = ode45(odefcn, tvec, x0);

tvec also provided a way of interpolating the model input vector inside the ODE:

ix = interp1(tvec, ix, t); % this is interpolating ix at t within tvec

After encountering difficulties with the discontinuities in the model, I am re-writing it to use event location so that I can stop and re-start the integration after making changes to the equations:

options = odeset('Events',@(t, x) EBDevents(t, x, ix, tvec, lr, krx, crx, adr, ...
    adr_dbz, Fxr_max));

[T, x, te, xe, ie] = ode45(odefcn, tvec, x0, options);

Eventually, I would like to end up with a concatenated set of integration results that align to tvec. So far it is working well for the first part of integration and I can now use the last returned values as initial conditions for the next part. However, I'm unsure what to do with tvec for the following periods of integration. The values returned in T and x are all at my requested intermediate points except the last one, which is at the event location. I guess I have the following options:

1. Use the last but one of the returned time and values as the start point for the next integration, plus the remaining points in tvec. I assume this would risk running into the same event the stopped the first integration.
2. Use the event location time and values (te and xe) as the start point and the points thereafter in tvec. After concatenating the results, there would be odd points that didn't align to tvec so some interpolation would be required.
3. Use just start and end times and no intermediate points for all the integration events and then interpolate with tvec after concatenating all the results. This wouldn't have exact points for the events but I can live with that. I could continue to pass tvec to the ODE for interpolating the input.

What's the most appropriate way forward?

Thanks, Simon.

Answer 1

[Jan]

The most accurate numerical method is to let the integrator decide for the best step sizes. Then you get the optimal solution between the accuracy through small step sizes and reducing the accumulated rounding errors by a small number of steps.

If you need the output of the trajectory to certain time points an interpolation after the integration is a good strategy. I did not take a look into the code of ODE45 for a long time now, but many other professional integrators do not modify the step size of the integration to meet requested time points for the output, but they re-use the calculated trajectory and the intermediate results from the Runge-Kutta steps for an interpolation.

This means, that your point 3. ist the best choice. To increase the accuracy, reduce the tolerances.

Another approach would be define a tspan to obatin the trajectories to defined time points, and to remove the event times afterwards:

% !!!UNTESTED CODE!!!
y0     = 0;
tspan0 = 1:1:10;
t0     = 1;
ready  = false;
Opt    = odeset('event', ???);
tAll   = [];
yAll   = [];
while ~ready
    tspan  = [t0, tspan0(tspan0 > t0)];  % [TYPO fixed]
    [t, y] = ode45(@fcn, tspan, y0, Opt);
    tAll   = cat(1, tAll, t);
    yAll   = cat(1, yAll, y);
    
    ready = (t(end) == tspan0(end));
    if ~ready
        y0 = y(end, :);
        t0 = t(end);
    end
end

Jitter          = ~ismembertol(tAll, tspan0);
tAll(Jitter)    = [];
yAll(Jitter, :) = [];

So crop tspan in each iteration and remove the irregular time points introduced by the event functions afterwards.

[Simon Aldworth]

Thanks Jan.

For the record, I think tspan should be:

tspan  = [t0, tspan0(tspan0 > t0)];

So that the last value of t returned from the most recent integration is concatenated with all the values in tspan0 higher than it.

Hope you agree.

Regards.

[Jan]

You are right, I've fixed this typo.