ABM stepper can't solve x'(t) = t accurate enough.

[GregorDeCillia]

If i try to solve this system with the adams_bashforth_multon solver of any order, I get errors around 0.01 when integrating from t = 0 to t = 1. That should not happen since the actual solution x(t) = t^2/2 is a polynomial, and the ABM stepper uses polynomial interpolations for integration.

Method used: integrate_const( abmStepper , rhs , x , 0.0 , 1.0 , 0.1 ) with void rhs(x,dx,t) {dx=t} Boost version: 1.57.0

I think the problem is the corrector method since the adams_bashforth stepper solves the exact same problem with accuracy 10^-16

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

typedef boost::array< double , 1 > state_type;
using namespace boost::numeric::odeint;

void f( state_type x, state_type &dxdt, double t ){
    dxdt[0] = t;
}

int main() {
    int nSteps = 10;
    double t;
    state_type x;

    t = 0.0;    x[0] = 0.0;
    integrate_const( adams_bashforth_moulton< 6 , state_type >(), 
                 f, x, t, 1.0, 1.0/nSteps );
    std::cout << x[0] << std::endl;      // 0.483507

    t = 0.0;    x[0] = 0.0;
    integrate_const( adams_bashforth< 6 , state_type >(), 
                f, x, t, 1.0, 1.0/nSteps );
     std::cout << x[0] << std::endl;     // 0.5
}

[GregorDeCillia]

Update: I found out that the right hand side never gets evaluated for t = 1 which should happen at the last corrector step. To see this I used the system function

void f( state_type x, state_type &dxdt, double t ){
    dxdt[0] = t;
    std::cout << "RHS evaluated ad t = " << t << std::endl;
}

Can anyone verify this? Btw, thanks for the great software. The runtime performance is really amazing.

[mariomulansky]

Thanks for your feedback! What happens if you use integrate_n_steps instead of integrate_const, your problem might be that integrate_const already stops at t=0.9 due to finite precision. Otherwise we have to check the abm for correctness.

[GregorDeCillia]

Thanks for the quick response.

Using the integrate_n_steps instead leads to the exact same results. Same if i try to make a for loop using sterpper.do_step. I am very sure the corrector evaluates the right hand side at the wrong time: t instead of t+dt

Here is what i just used for testing

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

typedef boost::array< double , 1 > state_type;
using namespace boost::numeric::odeint;

void f( state_type x, state_type &dxdt, double t ){
    dxdt[0] = t;
    std::cout << t << " ";
}

int main() {
    int nSteps = 10;
    double t;
    state_type x;

    t = 0.0;    x[0] = 0.0;    
    std::cout << "ABM\n\tRHS evaluations at t = " ;
    integrate_n_steps( adams_bashforth_moulton< 6 , state_type >(), 
               f, x, t, 1.0/nSteps, nSteps );
    std::cout  <<  "\n\tx(t=1) =" <<  x[0] << std::endl;      // 0.483507

    t = 0.0;    x[0] = 0.0;
    std::cout << "AB\n\tRHS evaluations at t = " ;
    integrate_const( adams_bashforth< 6 , state_type >(), 
                f, x, t, 1.0, 1.0/nSteps );
    std::cout << "\n\tx(t=1) = " << x[0] << std::endl;     // 0.5
}

[mariomulansky]

You were right, the corrector step was called with the wrong time value. This was not caught by our tests as we only tested with time-independent ODEs. I've fixed it and added a test case. Thanks a lot for recognizing this ugly bug and reporting here. Let me know if there are any more problems!

[GregorDeCillia]

Thanks!