How to programatically stop ode solver execution after fixed amount of computation time?

I have an initial value problem whose parameters I want to tune to measurement data. I found that my problem seems to become stiff if one or more parameters have very high values, at least what I observe is that in those cases, the ode45 solver takes a very long time and in most cases I end up terminating the calculation before it finishes. I tried to use instead the ode15s solver, which came to a solution in very short time.

I would like to automate this process somehow, and the easiest solution I can think of is to run the ode45 solver and after a given amount of time, terminate the execution and make a second attempt using the ode15s solver. I am aware of the existence of event triggers, but never worked with them and can't really figure out how I could use them in the context of my problem.

Can someone please draft how to achieve this? Thank you very much!

[Torsten]

I suggest using ODE15S throughout.

[Niklas Barthel]

Hey, thanks for your answer.

Could you point out what would be the drawbacks of doing as you suggest in cases where the parameters have values in the "non-stiff" regime, ie. when the ode45 solver would actually be able to solve the problem in a reasonable amount of time?

[Torsten]

No drawbacks concerning precision.

For higher-dimensional problems which are not stiff, using a stiff solver might increase the runtime and the required RAM because of the Jacobian matrix that has to be built in each integration step.

But changing the solver during integration is nonsense in my opinion.

[Niklas Barthel]

Thanks for pointing this out. I have accepted @Jan|s answer because it precisely answers my question, but I think I will in fact adapt your suggestion and use ode15s instead. Again, thanks!

Answer 1

[Jan]

The Event function is a bad idea: If the specified time has come, ode45 tries to find the simulated t, when the event value has changed. This wastes computing time.

The OutputFcn is thought for triggering an external event:

stopTime = now + 2 / 86400;  % 2 seconds
options  = odeset('OutputFcn', @(t, y, flag) myOutputFcn(t, y, flag, stopTime));
Sol      = ode45(@fcn, [0, 10], 0, options);

Sol.x(end) < 10   % Stopped by timer?

function dy = fcn(t, y)
dy = sin(t);
v  = exp(exp(sqrt(rand(1000))))^3;  % Waste time to slow down processing
end

function status = myOutputFcn(t, y, flag, stopTime)
status = double(now > stopTime);
end

This stops the integration after 2 seconds. This a kind of crash driven programming: instead of analysing the stiffness of the function to be intergrated, the runtime of the integrator is used to determine problems. Running this on a very slow or very fast machine can influence the method, so this is a bad idea. E.g. if Windows decides to install an update in the background, your algorithm might detect a wrong stiffness limit. Then repeated computations can cause different results.

[Niklas Barthel]

Hey, thanks for your answer!

"This a kind of crash driven programming: instead of analysing the stiffness of the function to be intergrated, the runtime of the integrator is used to determine problems." - This is something I am completely aware of, anyway thanks for pointing that out. This is meant to be a quick and dirty solution for the moment, and considering that the problem is very high-dimensional and thus I believe I would have a hard time analysing the stiffness of the problem, plus this is code that is - at least for now - only run on my own machine, I think this might be suitable for the moment.

I will accept your answer as soon as I implemented and briefly tested your code. Thanks again for your time.

[Jan]

The output function is called after an successfully accepted step. Then you can use the provided t,y to calculate the sensitivity matrix: The cumulative product of the matrix, which measures then quotient of the output of the function to be integrated to a variation of the input. In other words this measures the sensitivity of the solution to noise. This can be used to determine a bad choice of tolerances and you can detect stiffness also.

As long, as you know about the "dirtiness" and document this exhaustively, everything is fine.

[Niklas Barthel]

Hey, I was able to adapt your solution and everything works fine, thanks a lot for your time!