integrate_adaptive and integrate_times give different answers for negative step size

Question

I'm using the odeint library in Boost. When using the integrate_adaptive function, the results are as expected. However, when using integrate_times, the ODE is evaluated at very different times that are outside the range of integration. This is a problem for me because my ODE is not defined for some of the values that it is being evaluated at.

The code below demonstrates the issue. The x values for which the ODE is evaluated are printed to the screen.

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

struct observe
{
    std::vector<std::vector<std::complex<double> > > & y;
    std::vector<double>& x_ode;

    observe(std::vector<std::vector<std::complex<double> > > &p_y, std::vector<double> &p_x_ode) : y(p_y), x_ode(p_x_ode) { };

    void operator()(const std::vector<std::complex<double> > &y_temp, double x_temp)
    {
        y.push_back(y_temp);
        x_ode.push_back(x_temp);
    }
};

class Direct
{
    std::complex<double> alpha;
    std::complex<double> beta;
    std::complex<double> R;
    std::vector<std::vector<std::complex<double> > > H0_create(const double y);

public:
    Direct(std::complex<double> p_alpha, std::complex<double> p_beta, double p_R) : alpha(p_alpha), beta(p_beta), R(p_R) { }

    void operator() (const std::vector<std::complex<double> > &y, std::vector<std::complex<double> > &dydx, const double  x)
    {
        std::vector<std::vector<std::complex<double> > > H0 = H0_create(x);

        for(int ii = 0; ii < 6; ii++)
        {
            dydx[ii] = 0.0;
            for(int jj = 0; jj < 6; jj++)
            {
                dydx[ii] += H0[ii][jj]*y[jj];
            }
        }
    }
};


std::vector<std::vector<std::complex<double> > > Direct::H0_create(const double x)
{
    std::complex<double> i = std::complex<double>(0.0,1.0);
    std::cout << x << std::endl;
    double U = sin(x*3.14159/2.0);
    double Ux = cos(x*3.14159/2.0);
    std::complex<double> S = alpha*alpha + beta*beta + i*R*alpha*U;

    std::vector<std::vector<std::complex<double> > > H0(6);
    for(int ii = 0; ii < 6; ii++)
    {
        H0[ii] = std::vector<std::complex<double> >(6);
    }

    H0[0][1] = 1.0;
    H0[1][0] = S;
    H0[1][2] = R*Ux;
    H0[1][3] = i*alpha*R;
    H0[2][0] = -i*alpha;
    H0[2][4] = -i*beta;
    H0[3][1] = -i*alpha/R;
    H0[3][2] = -S/R;
    H0[3][5] = -i*beta/R;
    H0[4][5] = 1.0;
    H0[5][3] = i*beta*R;
    H0[5][4] = S;

    return H0;
}

int main()
{
    int N = 10;

    double x0 = 1.0;
    double xf = 0.0;

    std::vector<double> x_ode(N);
    double delta_x0 = (xf-x0)/(N-1.0);
    for(int ii = 0; ii < N; ii++)
    {
        x_ode[ii] = x0 + ii*delta_x0;
    }
    x_ode[N-1] = xf;

    std::vector<std::vector<std::complex<double> > > y_temp;
    std::vector<double> x_temp;

    std::complex<double> i = std::complex<double>(0.0,1.0);
    std::complex<double> alpha = 0.001*i;
    double beta = 0.45;
    double R = 500.0;
    std::complex<double> lambda = -sqrt(alpha*alpha + beta*beta + i*R*alpha);

    // Define Initial Conditions
    std::vector<std::complex<double> > ICs = {1, lambda, -i*alpha/lambda,0,0,0};

    // Initialize ODE class
    Direct direct(alpha,beta,R);

    {
        using namespace boost::numeric::odeint;

        double abs_tol = 1.0e-10;
        double rel_tol = 1.0e-6;

        std::cout << "integrate_adaptive x values:\n";
        size_t steps1 = integrate_adaptive(make_controlled<runge_kutta_cash_karp54<std::vector<std::complex<double> > > >(abs_tol, rel_tol), direct, ICs, x0, xf, delta_x0, observe(y_temp,x_temp));

        std::cout << "\n\nintegrate_times x values:\n";
        size_t steps2 = integrate_times(make_controlled<runge_kutta_cash_karp54<std::vector<std::complex<double> > > >(abs_tol, rel_tol), direct, ICs, x_ode.begin(), x_ode.end(), delta_x0, observe(y_temp,x_temp));
    }

    return 0;
}

I am compiling and running by using these commands:

g++ main.cpp -std=C++11; ./a.out

The code produces this output:

integrate_adaptive x values:
1
0.977778
0.966667
0.933333
0.888889
0.902778
0.888889
0.849758
0.830193
0.771496
0.693235
0.717692
0.693235
0.654104
0.634539
0.575842
0.497581
0.522037
0.497581
0.45845
0.438885
0.380188
0.301927
0.326383
0.301927
0.262796
0.24323
0.184534
0.106273
0.130729
0.106273
0.0850181
0.0743908
0.042509
0
0.0132841


integrate_times x values:
1
0.977778
0.966667
0.933333
0.888889
0.902778
0.888889
0.84944
0.829716
0.770543
0.691645
0.716301
0.777778
0.738329
0.718605
0.659432
0.580534
0.60519
0.666667
0.627218
0.607494
0.54832
0.469423
0.494078
0.555556
0.512422
0.490855
0.426154
0.339886
0.366845
0.444444
0.397898
0.374625
0.304806
0.211714
0.240805
0.333333
0.281908
0.256196
0.179058
0.0762077
0.108348
0.222222
0.170797
0.145085
0.0679468
-0.0349035
-0.00276275
0.111111
0.059686
0.0339734
-0.0431643
-0.146015
-0.113874
0.111111
0.0671073
0.0451054
-0.0209003
-0.108908
-0.0814056

The range of integration is from x = 1 to 0 but the ODE is being evaluated at x values less than 0 when using integrate_times.

Hmm, it looks like a bug in integrate_times. The solution is observed at all steps the stepper does not only at the time points defined by x_ode. I have to check it once more, but it should only "observe" the ode at 10 time points. Thank you for reporting this. — headmyshoulder, Oct 01 '13 at 20:44
@headmyshoulder Thanks for the quick reply. I really like your library because it is the best ODE solver that I've found for C++. For some reason I only see the problem when the state type is complex. Different systems that I'm working with that have state type of double work as expected. This isn't a huge problem for me though. Right now I am using the integrate_adaptive version and interpolating the solution onto the grid that I need. — OSE, Oct 01 '13 at 20:51
@headmyshoulder: the output happens in the rhs function, not in the observer. So the problem is only that the integrate_times overshoots the end point. This could be a bug related to going "backwards", i.e. negative dt. I will look into this — mariomulansky, Oct 01 '13 at 23:18
Ahh, you are right. I think you can add this point to your answer :) — headmyshoulder, Oct 02 '13 at 05:57
right, also it is much more efficient to use a dense-output stepper for integrate_times - I'll edit my answer. — mariomulansky, Oct 02 '13 at 15:44

mariomulansky · Accepted Answer · 2013-10-02T15:47:52.220

4

This is a bug in odeint due to the negative timesteps in your problem, I have created an issue on github: https://github.com/headmyshoulder/odeint-v2/issues/99 and I have implemented a fix. Please check out the latest odeint version from github and see if the problem remains. If so - feel free to open a new issue on github.

Thanks for pointing out that problem - and sorry for the bug.

Another note: I would suggest to use a dense-output stepper for the integrate_times routine as this is much more efficient (factor 2 in the best case). It basically does what you implemented as a fix above: using adaptive time-steps and interpolates at the intermediate points as required.

edited Oct 02 '13 at 15:47

answered Oct 02 '13 at 01:03

mariomulansky

933
5
7

Thanks for fixing this so quickly! – OSE Oct 02 '13 at 15:27
no problem - please consider using a dense-output stepper as I describe in my edited answer! – mariomulansky Oct 02 '13 at 15:48
1

Thanks for the tip, using a dense-output stepper really speeds things up. I have since switched to using the Dormand-Prince 5 algorithm because I had trouble getting dense-output to work with the Cash-Karp method. – OSE Oct 03 '13 at 20:16

integrate_adaptive and integrate_times give different answers for negative step size

1 Answers1