Closed nblinov closed 10 years ago
Hi, I fixed some small bugs in rosenbrock4_controller.hpp and commited it into the master branch.
Furhtermore you need to change the line with integrate_const to
size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< value_type > >( 1.0e-6 , 1.0e-6 ) , make_pair( stiff_system() , stiff_system_jacobi() ) , x , static_cast< value_type >( 0.0 ) , static_cast< value_type
( 50.0 ) , static_cast< value_type >( 0.01 ) , cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
Can you check if it is working for you?
On 09.11.2013 04:51, nblinov wrote:
|#include
include
include
include "../odeint-v2/boost/numeric/odeint.hpp"
include <boost/multiprecision/cpp_dec_float.hpp>
include <boost/phoenix/core.hpp>
include <boost/phoenix/core.hpp>
include <boost/phoenix/operator.hpp>
using namespace std; using namespace boost::numeric::odeint; namespace phoenix = boost::phoenix;
typedef boost::multiprecision::cpp_dec_float_50 value_type; typedef boost::array< value_type , 3 > state_type; //typedef double value_type;
//[ stiff_system_definition
typedef boost::numeric::ublas::vector< value_type > vector_type; typedef boost::numeric::ublas::matrix< value_type > matrix_type;
struct stiff_system { void operator()( const vector_type &x , vector_type &dxdt , value_type /* t / ) { dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 \ x[ 1 ]; dxdt[ 1 ] = x[ 0 ]; } };
struct stiff_system_jacobi { void operator()( const vectortype & /* x / , matrix_type &J , const valuetype & / t */ , vector_type &dfdt ) { J( 0 , 0 ) = -101.0; J( 0 , 1 ) = -100.0; J( 1 , 0 ) = 1.0; J( 1 , 1 ) = 0.0; dfdt[0] = 0.0; dfdt[1] = 0.0; } }; //]
int main( int argc , char **argv ) { //[ integrate_stiff_system vector_type x( 2 , 1.0 );
size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< value_type > >( 1.0e-6 , 1.0e-6 ) , make_pair( stiff_system() , stiff_system_jacobi() ) , x , 0.0 , 50.0 , 0.01 , cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" ); //] clog << num_of_steps << endl; return 0;
}|
I've added a test case for rosenbrock4+multiprecision and it runs fine with 56b4613. However, with the Rosenbrock stepper itself it is not possible to gain further accuracy by using multiprecision types because the internal constants of the stepper are only given with double precision. So you will not get a solution better than double precision with this stepper! This can only be overcome if you change the stepper to use constants defined with higher precision, but it might be even hard to find the constants with that accuracy in the literature. For the dopri5 stepper, for example, this is different as the internal constants there are given as rationals and thus in principle exactly and one benefits from using higher precision types!
Karsten, thanks, the example now runs fine. I wonder why you need to explicitly type cast though, can't boost::multiprecision do this implicitly?
Mario, yeah thanks for pointing that out. I also noticed that the constants were just written down as floats. I only need the extra precision for stability of the solution - increasing the working precision seemed to help out the mathematica stiff solver a bit for a similar system, so this is why i'm trying to get this working. I guess I'll have to dig around for more precise coefficients...
Thanks again!
Hello,
I'm trying to combine a stiff solver with boost::multiprecision. Changing stiff_system.cpp example to
gives rise to the following (terribly long error)
Any idea as to what's going on? Thanks for your help.