integrate_adaptive和integrate_times对于负步长给出不同的答案

Question

我在 Boost 中使用 odeint 库。使用integrate_adaptive函数时，结果符合预期。但是，使用Integrate_times 时，ODE 会在非常不同的时间进行计算，这些时间超出了积分范围。这对我来说是一个问题，因为我的 ODE 没有针对某些正在评估的值进行定义。

下面的代码演示了这个问题。计算 ODE 的 x 值将打印到屏幕上。

#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;
}

我正在使用以下命令进行编译和运行：

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

该代码产生以下输出：

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

积分的范围是从 x = 1 到 0，但使用Integrate_times 时，ODE 是在 x 值小于 0 时计算的。

Answer 1

由于问题中的负时间步长，这是 odeint 中的一个错误，我在 github 上创建了一个问题：https: //github.com/headmyshoulder/odeint-v2/issues/99 并且我已经实现了修复。请从github上查看最新的odeint版本，看看问题是否依然存在。如果是这样 - 请随意在 github 上打开一个新问题。

感谢您指出这个问题 - 并对这个错误表示歉意。

另一个注意事项：我建议对Integrate_times例程使用密集输出步进器，因为这样效率更高（最好情况下为因子2）。它基本上执行您在上面作为修复实现的操作：使用自适应时间步长并根据需要在中间点进行插值。