Adam Bashforth Moulten Inconsistent Results

[jbmccready]

Hello,

The adam_bashforth_moulten stepper seems to not be working properly. Using the simple harmonic oscillator example I get incorrect and inconsistent results. The adam_moulten stepper also fails to compile and adam_bashforth exits at run time when given a step number greater than 6 (document suggests up to 8 steps is allowed: http://www.boost.org/doc/libs/1_55_0/libs/numeric/odeint/doc/html/boost/numeric/odeint/adams_bashforth.html).

I'm using Boost 1.55.0 library, compiling Linux gcc 64 bit, and I have read the stepper docs extensively at: http://www.boost.org/doc/libs/1_55_0/libs/numeric/odeint/doc/html/index.html

Here's the test code, it's just the simple harmonic oscillator example modified:

#include <iostream>
#include <vector>

#include <boost/numeric/odeint.hpp>

typedef std::vector< double > state_type;

const double gam = 0.15;

void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
{
    dxdt[0] = x[1];
    dxdt[1] = -x[0] - gam*x[1];
}

int main(int /* argc */ , char** /* argv */ )
{
    using namespace std;
    using namespace boost::numeric::odeint;

    runge_kutta_dopri5 < state_type > dopri5;
    runge_kutta4 < state_type > rk4;
    adams_bashforth_moulton < 6 , state_type > abm6;
    adams_bashforth_moulton < 5 , state_type > abm5;
    adams_bashforth_moulton < 3 , state_type > abm3;
    adams_bashforth < 6 , state_type > ab6;
    adams_bashforth < 4 , state_type > ab4;
    adams_bashforth < 3 , state_type > ab3;

    adams_bashforth < 8 , state_type > ab8; // 8 steps, exits at run time (boost assert about algebra) see bottom
    adams_moulton < 4 , state_type > am4; // fails to compile (cannot find do_step method) see bottom

    state_type x(2);
    x[0] = 1.0;
    x[1] = 0.0;

    vector<state_type> x_vec;
    vector<double> times;

    integrate_const( dopri5 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "dopri5: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( rk4 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "rk4: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( abm6 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "abm6: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( abm6 , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
    cout << "abm6 more steps: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    double time = 0.0;
    abm6.initialize(harmonic_oscillator , x , time , 0.1);
    while ( time < 10.0) {
        abm6.do_step(harmonic_oscillator , x , time , 0.1);
        time += 0.1;
    }
    cout << "abm6 init and loop: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    time = 0.0;
    x_vec.clear();
    times.clear();
    while ( time < 10.0) {
        abm6.do_step(harmonic_oscillator , x , time , 0.1);
        time += 0.1;
    }
    cout << "abm6 loop: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( abm5 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "abm5: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( abm3 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "abm3: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( ab6 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "ab6: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( ab4 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "ab4: " << x[0] << ' ' << x[1] << endl;

    x[0] = 1.0; // reset
    x[1] = 0.0;
    integrate_const( ab3 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
    cout << "ab3: " << x[0] << ' ' << x[1] << endl;

      //this causes boost assert and exits
//    x[0] = 1.0; // reset
//    x[1] = 0.0;
//    integrate_const( ab8 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
//    cout << "ab8: " << x[0] << ' ' << x[1] << endl;

    //this fails to compile
//    x[0] = 1.0; // reset
//    x[1] = 0.0;
//    integrate_const( am4 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
//    cout << "am4: " << x[0] << ' ' << x[1] << endl;
}

The output for me is:

dopri5: -0.421909 0.246408
rk4: -0.421912 0.246405
abm6: 0.149653 -0.291708
abm6 more steps: 0.119897 -0.206122
abm6 init and loop: 0.229627 -0.233711
abm6 loop: 0.0910368 -0.330448
abm5: 0.154713 -0.295906
abm3: 0.159841 -0.304814
ab6: -0.421909 0.246409
ab4: -0.422013 0.246291
ab3: -0.420533 0.245226

I do not understand the ABM results, I have a much more complicated ODE system that has inconsistencies using ABM as well. Am I misusing the stepper? After reading the stepper page docs (http://www.boost.org/doc/libs/1_55_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html) I would have thought all three methods (integrate function, init and loop do_step, and straight loop do_step where do_step automatically inits on first few steps) should work the same. All three methods are giving incorrect results.

Thanks for the great library! I'm completely new to git hub, let me know if you need any more information, I am also interested in possibly implementing an adaptive error controlled ABM stepper (the ODE system I have is RHS computationally intensive and I'm doing millions of integrations).

-Buck

[headmyshoulder]

Hi,

thank you for reporting this issue. It seems that the initialization of the ABM mehtod does not work as expected. I will have a detailed look at it and come back to you within the next few days.

On 12/08/2013 09:53 PM, jbmccready wrote:

Hello,

The adam_bashforth_moulten stepper seems to not be working properly. Using the simple harmonic oscillator example I get incorrect and inconsistent results. The adam_moulten stepper also fails to compile and adam_bashforth exits at run time when given a step number greater than 6 (document suggests up to 8 steps is allowed: http://www.boost.org/doc/libs/1_55_0/libs/numeric/odeint/doc/html/boost/numeric/odeint/adams_bashforth.html).

I'm using Boost 1.55.0 library, compiling Linux gcc 64 bit, and I have read the stepper docs extensively at: http://www.boost.org/doc/libs/1_55_0/libs/numeric/odeint/doc/html/index.html

Here's the test code, it's just the simple harmonic oscillator example modified:

include

include

include <boost/numeric/odeint.hpp>

typedef std::vector< double > state_type;

const double gam = 0.15;

void harmonic_oscillator( const state_type &x , statetype &dxdt , const double /* t / ) { dxdt[0] = x[1]; dxdt[1] = -x[0] - gam_x[1]; }

int main(int /* argc / , char* /* argv */ ) { using namespace std; using namespace boost::numeric::odeint;

runge_kutta_dopri5 < state_type > dopri5;
runge_kutta4 < state_type > rk4;
adams_bashforth_moulton < 6 , state_type > abm6;
adams_bashforth_moulton < 5 , state_type > abm5;
adams_bashforth_moulton < 3 , state_type > abm3;
adams_bashforth < 6 , state_type > ab6;
adams_bashforth < 4 , state_type > ab4;
adams_bashforth < 3 , state_type > ab3;

adams_bashforth < 8 , state_type > ab8; // 8 steps, exits at run time (boost assert about algebra) see bottom
adams_moulton < 4 , state_type > am4; // fails to compile (cannot find do_step method) see bottom

state_type x(2);
x[0] = 1.0;
x[1] = 0.0;

vector<state_type> x_vec;
vector<double> times;

integrate_const( dopri5 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "dopri5: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( rk4 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "rk4: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( abm6 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "abm6: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( abm6 , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
cout << "abm6 more steps: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
double time = 0.0;
abm6.initialize(harmonic_oscillator , x , time , 0.1);
while ( time < 10.0) {
    abm6.do_step(harmonic_oscillator , x , time , 0.1);
    time += 0.1;
}
cout << "abm6 init and loop: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
time = 0.0;
x_vec.clear();
times.clear();
while ( time < 10.0) {
    abm6.do_step(harmonic_oscillator , x , time , 0.1);
    time += 0.1;
}
cout << "abm6 loop: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( abm5 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "abm5: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( abm3 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "abm3: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( ab6 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "ab6: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( ab4 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "ab4: " << x[0] << ' ' << x[1] << endl;

x[0] = 1.0; // reset
x[1] = 0.0;
integrate_const( ab3 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );
cout << "ab3: " << x[0] << ' ' << x[1] << endl;

  //this causes boost assert and exits

// x[0] = 1.0; // reset // x[1] = 0.0; // integrate_const( ab8 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 ); // cout << "ab8: " << x[0] << ' ' << x[1] << endl;

//this fails to compile

// x[0] = 1.0; // reset // x[1] = 0.0; // integrate_const( am4 , harmonic_oscillator , x , 0.0 , 10.0 , 0.1 ); // cout << "am4: " << x[0] << ' ' << x[1] << endl; }

The output for me is:

dopri5: -0.421909 0.246408 rk4: -0.421912 0.246405 abm6: 0.149653 -0.291708 abm6 more steps: 0.119897 -0.206122 abm6 init and loop: 0.229627 -0.233711 abm6 loop: 0.0910368 -0.330448 abm5: 0.154713 -0.295906 abm3: 0.159841 -0.304814 ab6: -0.421909 0.246409 ab4: -0.422013 0.246291 ab3: -0.420533 0.245226

I do not understand the ABM results, I have a much more complicated ODE system that has inconsistencies using ABM as well. Am I misusing the stepper? After reading the stepper page docs (http://www.boost.org/doc/libs/1_55_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html) I would have thought all three methods (integrate function, init and loop do_step, and straight loop do_step where do_step automatically inits on first few steps) should work the same. All three methods are giving incorrect results.

Thanks for the great library! I'm completely new to git hub, let me know if you need any more information, I am also interested in possibly implementing an adaptive error controlled ABM stepper (the ODE system I have is RHS computationally intensive and I'm doing millions of integrations).

-Buck

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

[mariomulansky]

i'm on that right now, the adams-bashforth now allows orders up to 8, and i've also added a accuracy test for adams-bashforth. i will look into the abm and am steppers now.

[headmyshoulder]

I'm on that too :)

The problem is in adams moulton. maybe the wrong coefficients...

On 09.12.2013 12:26, Mario Mulansky wrote:

i'm on that right now, the adams-bashforth now allows orders up to 8, and i've also added a accuracy test for adams-bashforth. i will look into the abm and am steppers now.

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

[jbmccready]

Thanks for the quick reply. For the abm stepper what is the purpose of using a combination of explicit (bashforth) and implicit (moulton) without any adaptive step sizing?

[headmyshoulder]

On 12/11/2013 03:54 AM, jbmccready wrote:

Thanks for the quick reply. For the abm stepper what is the purpose of using a combination of explicit (bashforth) and implicit (moulton) without any adaptive step sizing?

It is a method, which evaluates the system function only once or twice but has a high order up to 6 or 8. Therefore, this kinds of steppers can be used for systems where the evaluation of the r.h.s is expensive. The explicit step is a predictor step and the implicit step is a corrector step. Step size control is possible but not as straightforward as for the Runge-Kutta like steppers.

I think I can come up this morning for a fix of the abm stepper.

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

[headmyshoulder]

Ok, I pushed a fix.

[mariomulansky]

Ok I've just pushed a fix activating the orders 7 and 8 for adams bashforth moulton. With this all your issues should be resolved. Please check your tests again and let us know if some problems persist - otherwise we would close the issue. thanks a lot for reporting this, the adams-bashforth-moulton was indeed embarassingly screwed up - sorry!

[jbmccready]

It is a method, which evaluates the system function only once or twice but has a high order up to 6 or 8. Therefore, this kinds of steppers can be used for systems where the evaluation of the r.h.s is expensive. The explicit step is a predictor step and the implicit step is a corrector step. Step size control is possible but not as straightforward as for the Runge-Kutta like steppers.

Ah right, I wasn't really thinking about predictor/corrector pair, I was thinking about the moulton method being used to compare against the bashford method to determine error for adaptive step sizing. To add adaptive step sizing to the abm method I think it would require reevaluating past steps using a single step method when a step goes outside of tolerance. A variable order cash karp method may be more effective even when the r.h.s. is very expensive. This is all kind of drifting off of the issue, but are you guys interested in adding more steppers or options (such as a variable order cash karp, adaptive step sizing for abm, dense output for more steppers, etc.) to the library? I'm interested in possibly contributing.

I'm somewhat busy with finals this week, but I'll test the fixed abm method on my system within the next few days. Thanks for the quick response and fix.

[mariomulansky]

Hi, we are always very happy about contributions. We already have step-size control, e.g. with Cash-Karp, Dopri5 and Rk78 steppers, for the dopri5 there is also dense output available. However, a nice addition would be an error control for the adams-bashforth-moulton. I think this can be implemented based on the current abm stepper by replacing the single corrector step performed now with an iteration of corrector steps until the desired precision is reached (the series of corrector steps should converge). I found a short description and exemplary implementation there: http://www.mymathlib.com/diffeq/adams/ Again, we would be more than happy about contributions!

Thanks, Mario

[headmyshoulder]

On 12/11/2013 06:53 PM, jbmccready wrote:

It is a method, which evaluates the system function only once or twice
but has a high order up to 6 or 8. Therefore, this kinds of steppers can
be used for systems where the evaluation of the r.h.s is expensive. The
explicit step is a predictor step and the implicit step is a corrector
step. Step size control is possible but not as straightforward as for
the Runge-Kutta like steppers.

Ah right, I wasn't really thinking about predictor/corrector pair, I was thinking about the moulton method being used to compare against the bashford method to determine error for adaptive step sizing. To add adaptive step sizing to the abm method I think it would require reevaluating past steps using a single step method when a step goes outside of tolerance.

Step size control with Adams Bashforth method usually works by rescaling the past steps. Therefore reevaluation is not needed. It should be simple to add to odeint.

A variable order cash karp method may be more effective even when the r.h.s. is very expensive. This is all kind of drifting off of the issue, but are you guys interested in adding more steppers or options (such as a variable order cash karp, adaptive step sizing for abm, dense output for more steppers, etc.) to the library? I'm interested in possibly contributing.

We would be very happy about any contributions, but especially for more steppers. If you need any help or guidance let us know.

I'm somewhat busy with finals this week, but I'll test the fixed abm method on my system within the next few days. Thanks for the quick response and fix.

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

[jbmccready]

Alright I tested it on my system using different orders 5, 6, 7, 8, all work as expected. Thanks for the quick fix!

Step size control with Adams Bashforth method usually works by rescaling the past steps. Therefore reevaluation is not needed. It should be simple to add to odeint.

What I meant by reevaluation was the rescaling of past steps, step - 1, step - 2, step - 3, etc. have to be recomputed using a single step method, at least that's what I have gathered in my research. Here's a pretty concise breakdown: http://www.math.usm.edu/lambers/mat461/lecture6.pdf

Cash and Karp's paper published in 1990 also explains this recomputation of past steps for adaptive multistep methods: http://www.elegio.it/mc2/rk/doc/p201-cash-karp.pdf (page 202)

I'm interested in possibly adding the Cash and Karp variable order method as outlined in that paper, as it may be very powerful for some systems with expensive r.h.s. that transition frequently between smooth and sharp regions.

The other stepper I had in mind was Gobacki and Shampine's 2-3 order runge kutta (http://en.wikipedia.org/wiki/Bogacki%E2%80%93Shampine_method). This method may be well suited to what I am working on which is computing millions to billions of initial conditions and then mapping discrete convergence points, similar to what is shown here: http://simple.wikipedia.org/wiki/Magnetic_pendulum. I am only interested in knowing the final point converged to and therefore a higher error is possibly tolerable (I'm not sure without testing it, as much of these numerical methods really require testing to be sure of their effectiveness given a particular system).

I may implement these on my own first, or I could just start looking at the code and structure of odeint. Thanks again for the very quick response and fix.

[mariomulansky]

this sounds really great and would be a significant improvement of odeint! I think the best way would be to first do a straight-forward implementation, and then we can help putting it in the odeint framework. However, it would be helpful if you already try to modularize it, e.g. separate the step-size control from the stepper.

Anyhow, I will close this issue as I think the original problems with ab/abm methods are resolved?