ODE 45 integration control

I've been having some issues with running ODE45 on a model that is very flat unless forced with a very short duration square wave, in which case it forms a spike and returns to the baseline. The adaptive time step of ODE45 is causing the time step to be larger than the interval the model is forced for, and skipping the spikes (unless by chance, one of the t values lies in the forcing interval).

Is there any way to control the internal time steps used by the solver? I would like to be able to supply a list of time values and ensure that ode45 evaluates my function at those points. Elsewhere it can do what it likes I just need to make sure the function is evaluated within the forcing interval to trigger the shorter step size.

Many thanks,

Harry

Answer 1

[Jan]

Matlab's ODE solves cannot handle discontinuities. So controlling the time steps will not help reliably. See: http://www.mathworks.com/matlabcentral/answers/59582#answer_72047

The only way to treat discontinuities correctly is to stop the integration and restart it.

[Harry]

Thanks a lot for your comment, and stopping and restarting is working fine. The thing is if by chance one of the internal calculation values happens to lie within the interval its forced for, then it spikes perfectly, I just want that to happen every time. It's the right hand side of the ODE which is discontinuous, the solution is continuous but nonsmooth. A max step size of 0.1 also works perfectly, it's just very slow (computation time is very important in this application) and I only need that step size at certain values.

[Jan]

I repeat, that Matlab's ODE integrators cannot handle right hand sides with discontinuities correctly. This is a technical limitation.

[Harry]

It's not supposed to handle discontinuities correctly but in many cases it can. This is one of them if I can get one of the time steps to lie within the forcing period. I already know this will fix the problem. The question was how to implement this fix.

[Jan]

I do not agree that the integrators "can" handle discontinuities. It is true, that you get sometimes a final value. But from the view point of a scientist working in the field of numerical computations I do not call this final value a "result". The numerical integration is a kind of measurement process based on a large number of steps, which cause round off errors and local discretization errors. As in all measurements, the final value is a "result" only, if the measurement error can be estimated and given e.g. by a standard deviation. For an integration this means, that a sensitivity analysis is required. Unfortunately Matlab's ODE integrators do not provide the accumulated Wronski-matrix, but you at least vary the initial conditions to find out, how sensitive the final value reacts to this.

Imagine the integration of the equations of motion for a pen staying perpendicular on its tip. Would you think, that the calculated trajectory reflects the reality in any kind? An analysis of the sensitivity would reveal the instability of the system directly.

Now take into accound, that a discontinuos right hand side let the matrix of sensitivity explode and the integration even stops or fails, when such a discontinuity occurres exactly during the function evaluations during a step. Therefore I insist, that you can get a final value with some luck, but you'd run the integrator outside its specifications and therefore there is no scientific evidence, that the results are meaningful.

The handling of discontinuities in the right hand side for integrations has been one of the topic of my diploma thesis. Unfortuantely it was written in German, so it is not of general use. But as long as you found out, that stopping and restarting works fine, this is a nearly perfect solution already. For the equations of motion of e.g. a multibody system an additional step is required: The change of the RHS causes a jump in the velocities, e.g. when a foot of a running robot touches the floor. This concerns even the other bodies, which are connected by joints. Letting the integrator running stupidly over such discontinuities works sometimes, but it is numerically instable and the results are junk.

So my conclusion: I cannot really help to solve your question. You can modify the code of ODE45 to solve your needs, but I would recommend not to do it due to the mentioned problems.

[Harry]

Sorry for the delayed response, for some reason MATLAB never told me I had an additional comment. I think the main difference here is that the severity of the discontinuity is quite small, I can approximate it with a continuous function that's still quite smooth (like a Gaussian instead of a square pulse), it's just where it's flat for so long the step size becomes very large and jumps over the small period in which the derivative isn't zero. If it's very flat between 500 and 1000, with a Gaussian function between 1000 and 1001, but the internal time steps 999 and 1002 satisfy the error tolerances, the solver totally misses the Gaussian function, despite there being no discontinuity. Between 500 and 999 however, the step size can be as big as it wants because nothing happens. The point here is that the same problem persists even if I run within the specifications.

Regardless of that, stopping and restarting has been an effective method so far (as has just replacing the whole thing with fixed step size forward Euler). We know if our results are meaningful because they match our data, but I get your point that it leaves our work open to criticism. Considering our spike instigates a travelling wave, each point in that wave is technically dependant on a travelling discontinuity so stopping and restarting at the relevant times becomes much more difficult (although I do get a 'solution' only restarting at the initiating spike).

Thanks again for your help and sorry this response came so late, I've only just been notified of your additional comment.

Answer 2

[Mischa Kim]

Harry, in

 ...
 tspan   = linspace(0,1,11);
 options = odeset('AbsTol',1e-8);
 [t,X]   = ode45(@ODE,tspan,IC,options);

for which the solver returns the solution at specific data points defined by tspan. You can further use the options vector with odeset to control other solver options such as tolerances, step size, etc.

[Harry]

Thanks for your help. The issue is that it doesn't return the solution at that point, it interpolates the solution at that point from its internal t values, which are chosen to fit the error tolerance. It's these internal values that I need greater control over rather than the output values. If by chance, one of these internal values lies within the period it's forced for then it works fine.