2

I want to solve a stiff ode with odeint. I was following this (rosenbrock4_dense_output stepper) but, my function can grow pretty quickly so I want to stop the integration if x(t)>xMAX.

In this question, they have a solution for it but since I'm new with c++ I don't know how to implement this when using a rosenbrock4_dense_output stepper.

I would like to see how to write this specifically for rosenbrock4_dense_output stepper.

Community
  • 1
  • 1
marRrR
  • 35
  • 4
  • Just last week I added an example that finds the exact points where some threshold is crossed: https://github.com/headmyshoulder/odeint-v2/blob/master/examples/find_crossing.cpp, It uses the iterator interface of odeint together with a find_if to stop when the threshold is crossed. Then it does a bisection to approximate the crossing point, but you might not need that. To that adapt for rosenbrock you just have to change the stepper and the rhs function. – mariomulansky Oct 25 '15 at 23:19
  • This seems to be doing exactly what I need, but I don't know how to make it work with the rosenbrock4_dense_output stepper. Just changing runge_kutta_dopri5 to rosenbrock4_dense_output and instead of state_type using vector_type: where vector_Type is defined: typedef boost::numeric::ublas::vector< double > vector_type; is not all that is required. Would you mind explaining how this would change for my specific case. – marRrR Oct 26 '15 at 13:18
  • You also need to adjust the rhs, e.g. auto ode = make_pair( stiff_system() , stiff_system_jacobi() ); - see also headmyshoulder's answer below – mariomulansky Oct 26 '15 at 15:41

1 Answers1

2

Currently, this is not easily possible with odeint. If you can use the range library here you can combine a for_each and a find_if algorithm.

Otherwise you need to write the loop yourself, which is also not that difficult and should be similar to this:

auto stepper = make_dense_output< rosenbrock4< double > >( 1.0e-12 , 1.0e-12 );
auto ode = make_pair( stiff_system() , stiff_system_jacobi() );

double t = 0.0;
double const end_time = 50.0;
double const dt = 0.01;
vector_type x( 2 , 1.0 );
const double y_min = 1.0;

stepper.initialize( x , t , dt );
cout << t << " " << x[0] << " " << x[1] << endl; // or some other output
t += dt;
while( t < end_time )
{
    if( t > stepper.current_time() )
    {
        // perform a real step
        stepper.do_step( ode );
    }
    else  
    {
        // perform a dense output step
        stepper.calc_state( t , x );
        cout << t << " " << x[0] << " " << x[1] << endl; // or some other output
        t += dt;
    }
    if( x[1] < y_min ) // or some other condition
    {
        cout << "Bound reached." << endl;
        break;
    }
}
headmyshoulder
  • 6,240
  • 2
  • 20
  • 27
  • Do the odeint iterators not work with rosenbrock_dense_out? Or are you referring to the problem that with find_if we can not access the trajectory? – mariomulansky Oct 26 '15 at 02:02
  • The second problem: with `boost::range::find_if()` we can not access the trajectory and do the observation. – headmyshoulder Oct 26 '15 at 02:22
  • Would you mind writing explicitly how for the code will look with the rosenbrock4_dense_output stepper since I've been reading the documentation from odeint but I don't seem to get it right. – marRrR Oct 26 '15 at 11:46
  • Ok, done. `stiff_system` , `stiff_system_jacobi`are taken from the stiff system example you are referring to. – headmyshoulder Oct 26 '15 at 13:13