[Odeint] Stiff System of differential equations
Hi all, I am currently working with a very large set of differential equations; it is a disease model with 20160 classes depending on which vaccine has been admnistered, history of infection, and age. To solve this I have been using the runge_kutta_fehlber78 controlled stepper with the integrate_const function. There I need to set the absolute_error, relative_error and the step_size. My issue is that the results from the solver seems to depend highly on the step_size and whether or not an observer is provided. Since the solver is adaptive I feel like the results should not depend on step_size and they should definitely not depend on the presence of an observer. My question is why would the results be related to the step_size or the observer? How would it be possible to get the same results regardless of the step_size and the presence of an observer? Could this be related to the system being stiff and if so how would it be possible to solve it (The jacobian is too complicated to calculate because of the population dynamics)? Any help would be appreciated. Many thanks, Jules
Hi Jules, On 13.01.2016 21:01, jules wrote:
Hi all,
I am currently working with a very large set of differential equations; it is a disease model with 20160 classes depending on which vaccine has been admnistered, history of infection, and age. To solve this I have been using the runge_kutta_fehlber78 controlled stepper with the integrate_const function. There I need to set the absolute_error, relative_error and the step_size.
My issue is that the results from the solver seems to depend highly on the step_size and whether or not an observer is provided. Since the solver is adaptive I feel like the results should not depend on step_size and they should definitely not depend on the presence of an observer.
How do you use the solver? Runge Kutta 78 is an adaptive stepper, but you need to put it into the controlled_runge_kutta stepper. Furthermore, it is not a dense output solvers. This means that the solver need to reach to the time points where you want to observe your solution. Therefore, it needs to change the step size from the adaptive scheme just to reach the appropriate observation time. Maybe should try the dopri5 or bulirsh_stoer stepper, since they can use dense output and they should not be sensitive on the observation time and the time step.
My question is why would the results be related to the step_size or the observer? How would it be possible to get the same results regardless of the step_size and the presence of an observer? Could this be related to the system being stiff and if so how would it be possible to solve it (The jacobian is too complicated to calculate because of the population dynamics)? Any help would be appreciated.
Many thanks, Jules
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Hi Karsten, Thank you for your answer.
How do you use the solver? Runge Kutta 78 is an adaptive stepper, but you need to put it into the controlled_runge_kutta stepper.
This is how I used the runge_kutta_fehlberg78 stepper initally.
void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } typedef runge_kutta_fehlberg78< model_type , double , model_type , double , vector_space_algebra > stepper; double abs_err = .1; double rel_err = .00001; double step_size = .02; int steps = integrate_const( make_controlled<stepper>( abs_err , rel_err ) , model2, x, 0.0 , 10.0 , step_size, Printer_Observer);
Running it this way works but it takes a long time (which is fine and is to be expected). The results though seem to be depending on the abs_err and rel_err and stepsize very heavily. How do I know at what point the stepsize, abs_err and rel_err are small enough that I'll get the "right" answer? Also it gets stuck at t=0 if I set step_size=1.0.
Furthermore, it is not a dense output solvers. This means that the solver need to reach to the time points where you want to observe your solution. Therefore, it needs to change the step size from the adaptive scheme just to reach the appropriate observation time.
This makes sense but I don't feel like that should cause the solver to return an entirely different answer. I guess that is just the way of numeric solutions to ODEs though.
Maybe should try the dopri5 or bulirsh_stoer stepper, since they can use dense output and they should not be sensitive on the observation time and the time step.
I've since tried running the bulirsch_stoer stepper like this:
void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } double abs_err = .1; double rel_err = .00001; double step_size = .02; bulirsch_stoer_dense_out< model_type , double , model_type , double , vector_space_algebra > stepper {abs_err, rel_err}; int steps = integrate_const(stepper, model2, x, 0.0 , 10.0 , step_size, Printer_Observer);
The program compiles and runs but it seems like it gets stuck at t=0.02 (10+ minutes without moving anywhere). I've also tried with the runge_kutta_dopri5 with dense output:
void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } double abs_err = .1; double rel_err = .00001; double step_size = .02; int steps = integrate_const( make_dense_output( abs_err , rel_err , runge_kutta_dopri5< model_type , double , model_type, double , vector_space_algebra >()) , model2, x, 0.0 , 10.0 , step_size, Printer_Observer);
Here it also gets stuck at t=.02 and if I try with a step_size=1.0 it will get stuck at t=0. It seems like these steppers are also affected by the step_size and that they fail in the same way that the runge_kutta_fehlberg78 stepper did. Did these fail because they were improperly instantiated? Is the problem somehow too large or too stiff? Many thanks for your continued help, Jules
On 14.01.2016 22:48, Jules Tamagnan wrote:
Hi Karsten,
Thank you for your answer.
How do you use the solver? Runge Kutta 78 is an adaptive stepper, but you need to put it into the controlled_runge_kutta stepper.
This is how I used the runge_kutta_fehlberg78 stepper initally.
void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } typedef runge_kutta_fehlberg78< model_type , double , model_type , double , vector_space_algebra > stepper; double abs_err = .1; double rel_err = .00001; double step_size = .02; int steps = integrate_const( make_controlled<stepper>( abs_err , rel_err ) , model2, x, 0.0 , 10.0 , step_size, Printer_Observer);
Running it this way works but it takes a long time (which is fine and is to be expected). The results though seem to be depending on the abs_err and rel_err and stepsize very heavily. How do I know at what point the stepsize, abs_err and rel_err are small enough that I'll get the "right" answer? Also it gets stuck at t=0 if I set step_size=1.0.
In what way depend the results on the step size? In you example, the step size that you pass to the integrate function is in principle the observation time step. But it is also the initial time step which might be in appropriate for your use case. Can you try to solve the above example with integrate times: vector< double > obs_times = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 }; double step_size = 0.02; // now it is really the initial step size integrate_times( make_controlled< stepper >( abs_err , rel_err ) , model2 , x , step_size , times.begin() , times.end() , Print_observer );
Furthermore, it is not a dense output solvers. This means that the solver need to reach to the time points where you want to observe your solution. Therefore, it needs to change the step size from the adaptive scheme just to reach the appropriate observation time.
This makes sense but I don't feel like that should cause the solver to return an entirely different answer. I guess that is just the way of numeric solutions to ODEs though.
Yes, this is in principle correct. The problem with the above controlled stepper and integrate_const that they need to reach exactly the observation points. Hence, if you want to observe at say t = 0.5 the controlled stepper will break the time step adaption scheme to reach exactly at t = 0.5. This is not the case for the dense output steppers. Here intermediate values of your ode are interpolated (with the same accuracy as the solver).
Maybe should try the dopri5 or bulirsh_stoer stepper, since they can use dense output and they should not be sensitive on the observation time and the time step.
I've since tried running the bulirsch_stoer stepper like this:
void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } double abs_err = .1; double rel_err = .00001; double step_size = .02; bulirsch_stoer_dense_out< model_type , double , model_type , double , vector_space_algebra > stepper {abs_err, rel_err}; int steps = integrate_const(stepper, model2, x, 0.0 , 10.0 , step_size, Printer_Observer);
The program compiles and runs but it seems like it gets stuck at t=0.02 (10+ minutes without moving anywhere).
Hmm, this really looks like you have a problem with your ODE. You might be right, it could be too stiff. For analysis purposes you could also try to unroll the inner integration look to see what is happening here in detail: while( t < end_time ) { observ( x , t ); controlled_step_result res; size_t trials = 0; do { res = stepper.try_step( system , x , t , dt ); ++trials; // add some analysis to see the step size, e.g. // cout << t << " " << dt << endl; } while( ( res == fail ) && ( trials < max_attempts ) ); if( trials == max_attempts ) throw std::runtime_error( "Too many steps" ); } If you system is stiff you should also look at the Rosenbrock method. It is much better suited for stiff problems. But its needs to calculate the Jacobian.
I've also tried with the runge_kutta_dopri5 with dense output:
void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } double abs_err = .1; double rel_err = .00001; double step_size = .02; int steps = integrate_const( make_dense_output( abs_err , rel_err , runge_kutta_dopri5< model_type , double , model_type, double , vector_space_algebra >()) , model2, x, 0.0 , 10.0 , step_size, Printer_Observer);
Here it also gets stuck at t=.02 and if I try with a step_size=1.0 it will get stuck at t=0.
It seems like these steppers are also affected by the step_size and that they fail in the same way that the runge_kutta_fehlberg78 stepper did. Did these fail because they were improperly instantiated? Is the problem somehow too large or too stiff?
Many thanks for your continued help, Jules _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
In what way depend the results on the step size? In you example, the step size that you pass to the integrate function is in principle the observation time step. But it is also the initial time step which might be in appropriate for your use case.
After running a few more tests. It seems like the solution does not vary when step_size is changed but sometimes the solver will freeze, as I described earlier. For example, with my original code that uses the runge_kutta_fehlberg78 the solver gives the same answer when stepsize < 0.7 but when it it larger it will freeze and give no answer at all.
Can you try to solve the above example with integrate times:
vector< double > obs_times = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 }; double step_size = 0.02; // now it is really the initial step size integrate_times( make_controlled< stepper >( abs_err , rel_err ) , model2 , x , step_size , times.begin() , times.end() , Print_observer );
I used integrate_times as you described above but it also froze after the .1 step. After having waited for a very long time I eventually get this error message """ terminate called after throwing an instance of 'boost::exception_detail::clone_impl<boost::exception_detail::error_info_injector<std::overflow_error> >' what(): Integrate adaptive : Maximal number of iterations reached. A step size could not be found. """
If you system is stiff you should also look at the Rosenbrock method. It is much better suited for stiff problems. But its needs to calculate the Jacobian.
The biggest issue is that the Jacobian is too complicated to calculate, it depends on some population dynamics and higher order interpolation. Is there no other way to solve stiff systems. I've read about other methods such as CVODE and LSODA, but couldn't figure out how to get those to work.
vector< double > obs_times = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 }; double step_size = 0.02; // now it is really the initial step size integrate_times( make_controlled< stepper >( abs_err , rel_err ) , model2 , x , step_size , times.begin() , times.end() , Print_observer );
I used integrate_times as you described above but it also froze after the .1 step. After having waited for a very long time I eventually get this error message """ terminate called after throwing an instance of 'boost::exception_detail::clone_impl<boost::exception_detail::error_info_injector<std::overflow_error> >' what(): Integrate adaptive : Maximal number of iterations reached. A step size could not be found. """
Might it be possible that you reached something like a singular point? If the stepper does not find a step size it usually hints to some problems with the ODE. Can you show me roughly your equations?
If you system is stiff you should also look at the Rosenbrock method. It is much better suited for stiff problems. But its needs to calculate the Jacobian.
The biggest issue is that the Jacobian is too complicated to calculate, it depends on some population dynamics and higher order interpolation.
Is there no other way to solve stiff systems. I've read about other methods such as CVODE and LSODA, but couldn't figure out how to get those to work.
These methods exists but not within odeint and there are currently not reachable from odeint. You might have a look at here: https://computation.llnl.gov/casc/sundials/description/description.html
participants (3)
-
jules -
Jules Tamagnan -
Karsten Ahnert