передать параметр для повышения интегратора

Question

Я прочитал эти ответы 1 , 2 о передаче параметра в odeint (Boost).Я пытался здесь сделать тот же процесс.но код дает неправильный ответ.Вот моя попытка: Это упрощенный код ускоренного расчета ляпуновского показателя осциллятора Лоренца.

#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include "gram_schmidt.hpp"

using namespace std;
using namespace boost::numeric::odeint;

typedef vector<double> dim1;

const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
const size_t n = 3;
const size_t num_of_lyap = 3;
const size_t N = n + n * num_of_lyap;


// system with out perturbation
struct lorenz
{
    void operator()( const dim1 &x , dim1 &dxdt , double t ) const
    {
        dxdt[0] = sigma * ( x[1] - x[0] );
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -b * x[2] + x[0] * x[1];
    }
};

// system with perturbation
void lorenz_with_lyap( const dim1 &x , dim1 &dxdt , double t )
{
    lorenz()( x , dxdt , t );

    for( int l=0 ; l<num_of_lyap ; ++l )
    {
        const double *pert = &x[3 + l * 3]; 
        double *dpert = &dxdt[3 + l * 3]; 
        dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
        dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
        dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
    }
}

//------------------------------------------------------------------------

int main( int argc , char **argv )
{
    const double dt = 0.01;

    dim1 x(N); 
    x[0] = 10.; x[1] = 10.0; x[2] = 5.0; // initial condition
    dim1 lyap(3);

    runge_kutta4<dim1, double, dim1, double, range_algebra> rk4;
    // perform 10000 transient steps
    integrate_n_steps(rk4, lorenz(), x, 0.0, dt, 10000);

    fill( x.begin()+n , x.end() , 0.0 );
    for( size_t i=0 ; i<num_of_lyap ; ++i ) x[n+n*i+i] = 1.0;
    fill( lyap.begin() , lyap.end() , 0.0 );

    double t = 0.0;
    size_t count = 0;
    while( true )
    {

        t = integrate_n_steps( rk4 , lorenz_with_lyap , x , t , dt , 100 );
        gram_schmidt< num_of_lyap >( x , lyap , n );
        ++count;

        if( !(count % 1000) )
        {
            cout << t;
            for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
            cout << endl;
        }
    }
    return 0;
}

сначала я определил эти классы

class lorenz
{
public:
    void operator()...

и

class lorenz_with_lyap
{
public:
    void operator()(const dim1 &x, dim1 &dxdt, double t) const
    {
        dxdt[0] = sigma * (x[1] - x[0]);
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -b * x[2] + x[0] * x[1];
        for (int l = 0; l < num_of_lyap; ++l)
        {
            const double *pert = &x[3 + l * 3];
            double *dpert = &dxdt[3 + l * 3];
            dpert[0] = -sigma * pert[0] + sigma * pert[1];
            dpert[1] = (R - x[2]) * pert[0] - pert[1] - x[0] * pert[2];
            dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
        }
    }
};

тогда я могу немного изменить его, чтобы определить параметр:

class lorenz
{
    double sigma;
public:
    lorenz(double si) : sigma(si) { }
    void operator()...

, но определение класса lorenz_with_lyap меняет результат, и я не знаю, где я делаю что-то не так.

У вас есть идеи?(Теперь это домашнее задание, и я собираюсь проверить результат этого кода другим методом).Спасибо.

Answer 1

После некоторой борьбы, я наконец исправил это, выложите это здесь, если кто-то, вероятно, заинтересован:

class lorenz
{
public:
    void operator()(const dim1 &x, dim1 &dxdt, double t) const
    {
        dxdt[0] = sigma * (x[1] - x[0]);
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -b * x[2] + x[0] * x[1];
    }
};

class lorenz_with_lyap 
{
    const double sig;
public:
    lorenz_with_lyap(double si) : sig(si) { }
    void operator()(const dim1 &x, dim1 &dxdt, double t) const
    {
        dxdt[0] = sig * (x[1] - x[0]);
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -b * x[2] + x[0] * x[1];
        for (int l = 0; l < num_of_lyap; ++l)
        {
            const double *pert = &x[3 + l * 3];
            double *dpert = &dxdt[3 + l * 3];
            dpert[0] = -sig * pert[0] + 10.0 * pert[1];
            dpert[1] = (R - x[2]) * pert[0] - pert[1] - x[0] * pert[2];
            dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
        }
    }
};
//------------------------------------------------------------------------
// in main
        lorenz_with_lyap lo(10.0);
        t = integrate_n_steps(rk4, lo, x, t, dt, 100);
        ...