Skip to content


Subversion checkout URL

You can clone with
Download ZIP


ODE System exceptions #9

ph03 opened this Issue · 29 comments

8 participants


I just found your greate lib and I am just figuring out, if I can use it for my project. I'm doing flow visualization of dataset with a finite domain. The odeint lib could be used to integrate stream / path-lines in this domains. However, since the domains are only finite it will always be possible that a streamline integration leaves the domain. This is easily determined in a FlowField class that samples the field. However, I don't know how to communicate this event to the odeint integrators. I think the best way would be to integrate some kind of exception system to the odeint lib, for example for out of bounds sampling, critical point detection or simply for generic errors related to ODE system evaluation. Would it be hard to implement support for this
in a generic way or is the library not designed to handle this kind of ODE systems?

Greets Janick


Hi Janick,

yesterday we discussed a similar problem. I think this kind of application needs a special integrate function which handles the boundaries. This function might be implemented be the use of a general policy class capable of checking the bounds, for example

integrate_const( stepper , system , x , t_start , t_end , dt , boundary_handler );

boundary_handler is then a functor being called after each step, similar to the observers but which might change the state of the ODE. I think this is doable, although it will need some time for us to implement this.

Just two questions: What kind of steppers do you use? And what happens if a particle or streamline leaves the domain?


Hi! Right, a bounday_handler concept sounds like a very clean solution for this problem.
I haven't used odeint yet because of this current limitation. In our implementations most of the time a RungeKutta 45 scheme with adaptive stepsize control is more than sufficient (I think we have it from "Numerical Recipies", yours is called Dopri5 if I'm correct).
I think the best thing to happen when a streamline reaches a point outside of the boundary is to discard this point outside of the boundary (not even tell it to the observer) and instead redo a step from the previous point until the boundary is hit. I think depending on the used stepper there can be more specialized ways to solve this problem (actually we sample the polynomial curve that is implicitly generated by the integrator to obtain a point on the boundary), but I think a general viable way is to do some kind of binary search that reduces the stepsize of new points such that the boundary is hit "exactly" / in an epsilon distance. It would be also nice if the information of this "boundary hit" event will be retrievable, since the t_end value will not be reached. Particle movement outside the domain is undefined as there is no vector field data given.
However, I'm wondering if a boundary_handler concept is too restictive as there are perhaps more general exceptions possible (e.g. the observer notices some event and wants the integration to stop / redo the step with different stepsize because an error criteria is not met). I therefore think a general policy class implemented by the user that offers a high level of control of the integrator (maybe by a hierachy of different exceptions) and that communicates it to the user should be a good and general solution.

Thanks for your effort in this project! I would love to test any solution you come up with in our project, just tell me what to clone :)

Greets J


Hmm, after reading about your problem I think that other solutions might also fit well for you.

You might implement a step-size controller proxy which uses a controlled_runge_kutta stepper with a dopri5 stepper and implement your boundary checking logic within this proxy. So, the proxy calls your controlled stepper, performs the boundary checks and rejects or accepts the current step size and adjust the step size for the next step.

Another point you might take into account is dense output which interpolates your solution at arbitrary points between two steps. So, if a streamline leaves your domain in a step you can interpolate the value of streamline at the boundary. But you must be careful, since the current step might have been obtained from informations from the exterior of your range.



is there any news on allowing events to control the integration process? I am currently using the adaptive integrator for tracing rays and I would like to both stop the integration at the boundaries of my model and to catch other non-terminal events, like an intercept of an interface. The observer should allow me to trigger on or at least register non-terminal events, but there is no way as far as I can see to indicate that the integration should stop.

Thanks for the great software.

Update: I realized from a discussion on the mailing list[1] that I could probably do something with iterators.

Best regards, Gaute




sorry there are no news about conditional integrators. It is still on our list. Actually I first want to finalize the iterators before implementing the conditional integrate. The Iterators could also solve your problem, at least in a limited fashion. But they are currently redesgined, because they had some serious problems with memory management. I think the iterators will be finished within the next week. I will keep you informed here.


First of all, thanks for the excellent library.

Second, I'm also looking for a simple way to end integration when some function of the states undergoes a zero crossing. The other two ODE integrators that I am familiar with that have this functionality are Matlab (ode45 et al.) and the LSODAR function that is part of ODEPACK. There are lots of uses cases for this functionality and it would be great if this were part of odeint. I know a number of people who use Matlab solely for this events functionality but don't have the time/inclination/skills to write this functionality.

Are there any efforts at the moment to add this functionality?


Thanks for pointing me to those examples.

Is there an example which illustrates how to find (to some specified accuracy) the time at which the zero crossing occurs? Also, it would be nice to have an example which allows for setting up the function whose zero crossings you want to find and/or terminate integration on, though I guess the example which uses x[0] > 0.0 ? true : false could be expanded to have more sophisticated logic than just checking one of the states.

Being able to specify the direction of the zero crossing of the function of the states is handy. For example allowing g(x, t) (the function whose zeros you want to find during integration) to go from positive to negative and not trigger the event, but when it goes from negative to positive, have it trigger the event.


Also, if there is anything I can do to help, let me know. I'm still getting comfortable with the library, my knowledge of C++ template is fairly low, but if you have things to test out or just want to bounce ideas around, feel free to do so.


I achieved events/exceptions (in ray tracing) by re-implementing the adaptive integrator with the event/exception logic in between, it is then fairly easy to re-iterate on an event to reach a given accuracy for the event.


I will factor it out, but it is not very generic.


Check out:

When an event is detected, e.g.: I passed the stepper and state to a function which would then continue to step it (depending on the problem) to get the desired accuracy for the event. This is because, in my case, I could often quickly determine whether an event had happened between two steps with the overall stepper accuracy, but needed more time to determine exactly when (and at which state) the event happened in the range between the two steps.

When reading, bear in mind that this code has been relatively quickly cut out from a larger code base and some comments/stuff doesn't necessarily seem logical here. The most interesting part is Ray::solve () in ray.cpp.


I also made an implementation which required the use of a stepper with dense output. I then used bisection to find the zero-crossing. Here is an example based on the lorenz-example from the docs using to stop the first time the trajectory crosses x = 0:



Hi guys,

Hehe, I have already implemented a first version of integrate_conditional some month ago. It is already quite advanced but not completed yet. Have a look into the integrate_conditional branch. It provides an integrate_conditional function which takes an additional controller argument. The controller will take care about stepping, stopping condition, and what will happen if the condition is fulfilled.

Have a look at the There are two controllers: conditional_stop which simply stops if a give condition has been fulfilled. The second is approximate which halves the stepsize if a given condition is met until a specified precision is met.

I think all your use case can easily be mapped into specific controllers, especially the case with the bisection. Can you check if this model would work for you? I would then start to finish the implementation of conditional_integrate, create some often used controllers, and write documentation. Of course, any help is appreciated.

Best regards,



Hi Karsten,

Can this handle multiple events in the integrate_conditional, like if you're detecting zero crossings in two state variables?

Thanks for the great functionality on this, it makes switching from MATLAB way easier!



Can you send us an short example how you use integrate_conditional?

I actually didn't attempt that yet. I just started porting most of my simulations, so sorry I didn't succeed in using integrate_conditional!

Instead we decided to go for the iterator solution which already works in the master branch.

I guess I will have to implement it this way then. Do you have example code somewhere that demonstrates this? I found this example but I am confused how stop the integration. Can you elaborate on that?

Hi I found this useful example on how to use pure iterators on the master branch but I am currently stuck trying on a compiler error:

/usr/include/boost/iterator/iterator_facade.hpp:327:49: error: no matching function for call to ‘implicit_cast(std::pair&, const double&>*)’
/usr/include/boost/iterator/iterator_facade.hpp:327:49: note: candidate is:
/usr/include/boost/implicit_cast.hpp:18:10: note: template T boost::implicit_cast(typename boost::mpl::identity::type)

whenever I try to dereference the iterator. Do you know what might be wrong here?

Apologies if this is the wrong place to ask that.



Thanks again for the snippet!

I have just one more question. Is there a way to add an observer to make_adaptive_range ? I want to get all the steps in the integration as well.

Maybe I need to use the example in line 160?

Thanks again!


My apologies if this is already supported by the current implementation; it would be nice if the conditional integrator accepted input from the return value of the condition checker so that a distance or direction could be returned. If, for instance, I want to trace a ray untill a current interface the condition could more easily be achieved if the integrator knows whether it should go backwards or forwards, or how far it is from the interface. In MATLAB [0] this is implemented by using a return value (positive or negative) where 0 is condition met (with a given error, of course). In cases where distance or direction does not apply, returning 1 or 0 (condition met) would be sufficient. The integrator should also return from what direction the condition was met.



Yes, in #9 (comment) I wrote an integrator based on the integrate_adaptive code. But the controlling code was all included in the integrating function. This integrator detects if the event happens between two points, then uses the integrator to further approach the event. A similar approach could probably be generalized to solve my use case above for more general problems.


For a project I am working on I need to use an iterator-based solution in order to stop the integration when certain conditions are met.
It looks like the following code, when applied to our specific situation, will work.

      auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type > () );
      state_type x = {{ 10.0 , 10.0 , 10.0 }};
      auto iter = boost::find_if( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                             []( const std::pair< const state_type & , double > &x ) 
                                             {return ( x.first[0] < 0.0 ); }

      cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";

We are also interested in using the Bulirsch-Stoer method, but when I try to implement the Bulirsch-Stoer stepper as follows:

    bulirsch_stoer<state_type> stepper(1e-9, 0.0, 0.0, 0.0);        
    state_type x = {{ 10.0 , 10.0 , 10.0 }};

    auto iter = boost::find_if( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                 []( const std::pair< const state_type & , double > &x ) {
                                    return ( x.first[0] < 0.0 ); } );

    cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";

I get the following compilation error.

1>c:\program files (x86)\microsoft visual studio 10.0\vc\include\xutility(275): error C2582: 'operator =' function is unavailable in 'boost::numeric::odeint::adaptive_time_iterator'
1> with
1> [
1> Stepper=boost::numeric::odeint::bulirsch_stoer,
1> System=lorenz,
1> State=state_type,
1> StepperTag=boost::numeric::odeint::controlled_stepper_tag
1> ]

Could anyone help me implement the Bulirsch-Stoer stepper in make_adaptive_time_range?

Thanks for the help!

@mariomulansky mariomulansky referenced this issue from a commit
@mariomulansky mariomulansky +example for event detection, see #9
simple example for an event detection implementation based on dopri5
stepper and simple bisection.
Simple, straight forward implementation, but shouldnt be difficult to adapt
to other problems.

Hi everyone,

After some discussion I've revisited this problem and added a short example based on dopri5, find_if and a final bisection: 766d2cb
Quite straight-forward and shouldnt be too difficult to adapt. I hope this is a helpful example for some.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.