headmyshoulder / odeint-v2

odeint - solving ordinary differential equations in c++ v2
http://headmyshoulder.github.com/odeint-v2/
Other
342 stars 101 forks source link

long double support in odeint #102

Closed Crakem closed 10 years ago

Crakem commented 11 years ago

Hi, I'm using your great library for integrating a stiff (solar system) problem requiring long double precision (I have tested with octave on amd64 but it's slow) but can't compile My code (following your guide) typedef long double dbl; ... controlled_stepper_type controlled_stepper(default_error_checker< dbl >( AbsTol , RelTol , a_x , a_dxdt ) ); size_t steps = integrate_adaptive( controlled_stepper , ss9 , y , t0 , t_final, step ,record_state_and_time(x_vec, times));

error was: main.cpp: In function 'int main(int, char**)': main.cpp:136:112: error: no matching function for call to 'boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector > >::controlled_runge_kutta(boost::numeric::odeint::default_error_checker)' /usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:173:5: note: candidates are: boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::controlled_runge_kutta(const error_checker_type&, const stepper_type&) [with ErrorStepper = boost::numeric::odeint::runge_kutta_cash_karp54<std::vector >, ErrorChecker = boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, Resizer = boost::numeric::odeint::initially_resizer, error_checker_type = boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, stepper_type = boost::numeric::odeint::runge_kutta_cash_karp54<std::vector >] /usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:115:7: note: boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector > >::controlled_runge_kutta(const boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector > >&) In file included from /usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:25:0, from /usr/include/boost/numeric/odeint/integrate/integrate.hpp:24, from /usr/include/boost/numeric/odeint.hpp:64, from SolarSystem.cpp:4, from main.cpp:10: /usr/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp: In function 'size_t boost::numeric::odeint::detail::integrate_adaptive(Stepper, System, State&, Time&, Time, Time&, Observer, boost::numeric::odeint::controlled_stepper_tag) [with Stepper = boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector > >, System = solarsim, State = std::vector, Time = long double, Observer = record_state_and_time, size_t = long unsigned int]': /usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:44:61: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(Stepper, System, State&, Time, Time, Time, Observer) [with Stepper = boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector > >, System = solarsim, State = std::vector, Time = long double, Observer = record_state_and_time, size_t = long unsigned int]' main.cpp:137:126: instantiated from here /usr/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp:98:13: error: no matching function for call to 'boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector > >::try_step(solarsim&, std::vector&, long double&, long double&)'

headmyshoulder commented 11 years ago

Can you show us some more code? It is not possible for me to figure out what is happening here.

On 10/13/2013 02:16 AM, Crakem wrote:

Hi, I'm using your great library for integrating a stiff (solar system) problem requiring long double precision (I have tested with octave on amd64 but it's slow) but can't compile My code (following your guide) typedef long double dbl; ... controlled_stepper_type controlled_stepper(default_error_checker< dbl >( AbsTol , RelTol , a_x , a_dxdt ) ); size_t steps = integrate_adaptive( controlled_stepper , ss9 , y , t0 , t_final, step ,record_state_and_time(x_vec, times));

error was: main.cpp: In function 'int main(int, char**)': main.cpp:136:112: error: no matching function for call to 'boost::numeric::odeint::controlled_runge_kutta >

::controlled_runge_kutta(boost::numeric::odeint::default_error_checker)' /usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:173:5: note: candidates are: boost::numeric::odeint::controlled_runge_kutta::controlled_runge_kutta(const error_checker_type&, const stepper_type&) [with ErrorStepper = boost::numeric::odeint::runge_kutta_cash_karp54 >, ErrorChecker = boost::numeric::odeint::default_error_checker, Resizer = boost::numeric::odeint::initially_resizer, error_checker_type = boost::numeric::odeint::default_error_checker, stepper_type = boost::numeric::odeint::runge_kutta_cash_karp54 >] /usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:115:7: note: boost::numeric::odeint::controlled_runge_kutta > ::controlled_runge_kutta(const boost::numeric::odeint::controlled_runge_kutta > >&) In file included from /usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:25:0, from /usr/include/boost/numeric/odeint/integrate/integrate.hpp:24, from /usr/include/boost/numeric/odeint.hpp:64, from SolarSystem.cpp:4, from main.cpp:10: /usr/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp: In function 'size_t boost::numeric::odeint::detail::integrate_adaptive(Stepper, System, State&, Time&, Time, Time&, Observer, boost::numeric::odeint::controlled_stepper_tag) [with Stepper = boost::numeric::odeint::controlled_runge_kutta > >, System = solarsim, State = std::vector, Time = long double, Observer = record_state_and_time, size_t = long unsigned int]': /usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:44:61: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(Stepper, System, State&, Time, Time, Time, Observer) [with Stepper = boost::numeric::odeint::controlled_runge_kutta > >, System = solarsim, State = std::vector, Time = long double, Observer = record_state_and_time, size_t = long unsigned int]' main.cpp:137:126: instantiated from here /usr/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp:98:13: error: no matching function for call to 'boost::numeric::odeint::controlled_runge_kutta > >::try_step(solarsim&, std::vector&, long double&, long double&)'

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/102.

headmyshoulder commented 11 years ago

Can you change the line

typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;

to

typedef controlled_runge_kutta< error_stepper_type , default_error_checker< dbl > > controlled_stepper_type;

and test once more?

Crakem commented 11 years ago

Sure! I have prepared a test case on linux, main.cpp

#include <cstdlib>
#include "raytracer.cpp"
int main(int argc, char** argv) {
    using namespace std;
    using namespace boost::numeric::odeint;
    const dbl xpos0=ro+1e3;

    state_type x(4);
    x[0] = 0.0;
    x[1] = c/xpos0;
    x[2] = xpos0;
    x[3] = 0.0;
    double step,dist,RelTol,AbsTol;
    argc>1?step=std::atof(argv[1]):step=2e-2;       //paso
    argc>2?dist=std::atof(argv[2]):dist=2e12;        //distancia de integracion
    argc>3?RelTol=std::atof(argv[3]):RelTol=1e-18;        //tol relativa
    argc>4?AbsTol=std::atof(argv[4]):AbsTol=1e-7;        //tol absoluta

    light_ray rt(0.0);

    typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
    typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
    controlled_stepper_type controlled_stepper(default_error_checker< dbl >( AbsTol ,  RelTol , 1.0 , 1.0 ) );
    size_t steps = integrate_adaptive( controlled_stepper , rt , x , 0.0 , ((dbl)dist)*ro/c, step);

    std::cout << "Paso: " << step << ", Dist: " << dist << " (" << ((dbl) dist) * ro / c << " seg), aTol: " << AbsTol << " rTol: " << RelTol << std::endl;
    std::cout << "Evaluaciones: " << steps << std::endl;
    return (EXIT_SUCCESS);
}

and class: raytracer.cpp,

#include <iostream>
#include <vector>
#include <boost/numeric/odeint.hpp>

typedef long double dbl;
typedef std::vector<dbl> state_type;

const dbl G= 6.67259e-11;
const dbl GMs=1.32712440017987e20;
const dbl c=299792458.0;
const dbl ro=6.963e8;

class light_ray{
    dbl m_cm;
public:
    light_ray(dbl cm):m_cm(cm){}
    void operator()(const state_type &x , state_type &dxdt , const double /* t */)
    {
       dbl p=x[1],r=x[2],q=x[3];
       const dbl rr[]={r*std::cos(x[0]),r*std::sin(x[0])};

       dxdt[0]=p;
       dxdt[1]=-q*2.0*p/r;
       dxdt[2]=q;
       dxdt[3]=r*p*p-GMs/(r*r);
    }
};

Thanks for replies. Best regards, Enrique

Crakem commented 11 years ago

Hi, I've updated my code like this

//probar stiff
    typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
    typedef controlled_runge_kutta< error_stepper_type , default_error_checker< dbl > > controlled_stepper_type;
    //typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
    controlled_stepper_type controlled_stepper(default_error_checker< dbl >( AbsTol ,  RelTol , 1.0 , 1.0 ) );
    size_t steps = integrate_adaptive( controlled_stepper , rt , x , 0.0 , ((dbl)dist)*ro/c, step /*,record_state_and_time(x_vec, times)*/);
//    ]

but now error was, main.cpp:39:139: error: no matching function for call to 'integrate_adaptive(main(int, char**)::controlled_stepper_type&, light_ray&, state_type&, double, dbl, double&)'

headmyshoulder commented 11 years ago

Can you try to explictly cast the values for the times to long double, i.e.

size_t steps = integrate_adaptive( controlled_stepper , rt , x , static_cast< long double >( 0.0 ) , static_cast< long double >( ((dbl)dist)*ro/c ) , static_cast< long double >( step ) /*,record_state_and_time(x_vec, times)*/);

Does this work?

Crakem commented 11 years ago

Error now was, /usr/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp: In function 'size_t boost::numeric::odeint::detail::integrate_adaptive(Stepper, System, State&, Time&, Time, Time&, Observer, boost::numeric::odeint::controlled_stepper_tag) [with Stepper = boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector >, boost::numeric::odeint::default_error_checker >, System = light_ray, State = std::vector, Time = long double, Observer = boost::numeric::odeint::null_observer, size_t = long unsigned int]': In file included from /usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:25:0, /usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:44:61: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(Stepper, System, State&, Time, Time, Time, Observer) [with Stepper = boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector >, boost::numeric::odeint::default_error_checker >, System = light_ray, State = std::vector, Time = long double, Observer = boost::numeric::odeint::null_observer, size_t = long unsigned int]' /usr/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp:81:110: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(Stepper, System, State&, Time, Time, Time) [with Stepper = boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector >, boost::numeric::odeint::default_error_checker >, System = light_ray, State = std::vector, Time = long double, size_t = long unsigned int]' main.cpp:39:230: instantiated from here /usr/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp:98:13: error: no matching function for call to 'boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector >, boost::numeric::odeint::default_error_checker >::try_step(light_ray&, std::vector&, long double&, long double&)'

mariomulansky commented 11 years ago

There are a few minor bugs in your main.cpp. I've corrected them and with the following main.cpp it compiles for me with gcc 4.7.1:

#include <cstdlib>
#include "raytracer.cpp"

int main(int argc, char** argv) {
    using namespace std;
    using namespace boost::numeric::odeint;
    const dbl xpos0=ro+1e3;

    state_type x(4);
    x[0] = 0.0;
    x[1] = c/xpos0;
    x[2] = xpos0;
    x[3] = 0.0;
    dbl step,dist,RelTol,AbsTol;
    argc>1?step=std::atof(argv[1]):step=2e-2;       //paso
    argc>2?dist=std::atof(argv[2]):dist=2e12;        //distancia de integracion
    argc>3?RelTol=std::atof(argv[3]):RelTol=1e-18;        //tol relativa
    argc>4?AbsTol=std::atof(argv[4]):AbsTol=1e-7;        //tol absoluta

    light_ray rt(0.0);

    typedef runge_kutta_cash_karp54< state_type , dbl > error_stepper_type;

    size_t steps = integrate_adaptive( make_controlled<error_stepper_type>(AbsTol ,  RelTol) ,
                                       rt , x , 
                                       static_cast< dbl >(0.0) , 
                                       static_cast< dbl >(((dbl)dist)*ro/c), 
                                       static_cast< dbl >(step) );

    std::cout << "Paso: " << step << ", Dist: " << dist << " (" << ((dbl) dist) * ro / c << " seg), aTol: " << AbsTol << " rTol: " << RelTol << std::endl;
    std::cout << "Evaluaciones: " << steps << std::endl;
    return (EXIT_SUCCESS);
}

you need an extra template arguments for the stepper so it uses long double also for the parameters. Note how make_controlled makes it easier to get a controlled stepper. Finally, let me just say that your constants in the raytracer.hpp are initialized only with double accuracy, not long double... just to let you know :)

Crakem commented 11 years ago

oh, I'm sorry, now I'm getting,

/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:363:17: error: no matching function for call to 'max(double, boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<std::vector<long double, std::allocator >, long double, std::vector<long double, std::allocator >, long double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<long double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_tag>::value_type&)'

please, could you told me which version are you using? I'm using: boost-1.53.0 gcc-4.5.2

Thanks a lot for so quick replies :)

mariomulansky commented 11 years ago

ok, now that's a bug in the boost version of odeint which has been fixed on github. if you take the latest odeint version it should work.

Crakem commented 11 years ago

Thanks a lot for your assistance and your great, great, great software. I'll test it in next update. My gentoo box show me boost 1.54.0 is current one, maybe this change will be committed in 1.55.0? I would like to say your odeint is the fastest way for solving ode's :):):)

headmyshoulder commented 11 years ago

On 10/25/2013 08:45 PM, Crakem wrote:

Thanks a lot for your assistance and your great, great, great software. I'll test it in next update. My gentoo box show me boost 1.54.0 is current one, maybe this change will be committed in 1.55.0?

We will not get it into 1.55.0, since the release branch already closes on Monday and we need to run all the tests until Monday. But we include it in 1.56.

I would like to say your odeint is the fastest way for solving ode's :):):)

Thank you :)

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/102#issuecomment-27116369.

Crakem commented 11 years ago

Hi, I have found svn boots repository and I'm thinking about making a live ebuild for my distro; please, could you told me which revision I have to commit for testing your update? Thanks!!

mariomulansky commented 11 years ago

no also the svn boost repos dont contain the latest odeint version. you have to get odeint itself from here (github) and put it somewhere on your machen. if you then include the standalone odeint before boost in your built it will be taken and your code should compile.