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