В качестве первого примера мы возьмем самый простой случай и реализуем наш собственный вектор<my_vector>, который обеспечит контейнерный интерфейс.Укрепление. Диапазонработает вне коробки. Мы добавляем немного функциональности к нашему вектору, что позволяет выделить некоторую пропускную способность по умолчанию. Это полезно при использовании изменения размера, поскольку тогда изменение размера может быть гарантировано, чтобы не требовать нового распределения.

template< size_t MAX_N >
class my_vector
{
    typedef std::vector< double > vector;
public:
    typedef vector::iterator iterator;
    typedef vector::const_iterator const_iterator;
public:
    my_vector( const size_t N )
        : m_v( N )
    {
        m_v.reserve( MAX_N );
    }
    my_vector()
        : m_v()
    {
        m_v.reserve( MAX_N );
    }
// ... [ implement container interface ]

Единственное, что нужно сделать, кроме определения, это объявить мой вектор размером:

// define my_vector as resizeable
namespace boost { namespace numeric { namespace odeint {
template<size_t N>
struct is_resizeable< my_vector<N> >
{
    typedef boost::true_type type;
    static const bool value = type::value;
};
} } }

If we wouldn't specialize the is_resizeable template, the code would still compile but odeint would not adjust the size of temporary internal instances of my_vector and hence try to fill zero-sized vectors resulting in segmentation faults! The full example can be found in my_vector.cpp

Если ваш тип состояния работает сBoost.Range, но обрабатывает изменение размера по-разному, вам необходимо специализироваться на двух реализациях, используемых odeint для проверки размера состояния и изменения размера:

В качестве примера мы будем использовать<std::list>как тип состояния в одейте. Поскольку<std::list>не поддерживается<boost::size>, мы должны заменить реализацию того же размера и размера, чтобы получить список для работы с odeint. Следующий код показывает необходимые специализации шаблонов:

typedef std::list< double > state_type;
namespace boost { namespace numeric { namespace odeint {
template< >
struct is_resizeable< state_type >
{ // declare resizeability
    typedef boost::true_type type;
    const static bool value = type::value;
};
template< >
struct same_size_impl< state_type , state_type >
{ // define how to check size
    static bool same_size( const state_type &v1 ,
                           const state_type &v2 )
    {
        return v1.size() == v2.size();
    }
};
template< >
struct resize_impl< state_type , state_type >
{ // define how to resize
    static void resize( state_type &v1 ,
                        const state_type &v2 )
    {
        v1.resize( v2.size() );
    }
};
} } }

С помощью этих определений одеинт знает, как изменять размер<std::list>с, и поэтому их можно использовать в качестве типов состояний. Полный пример можно найти вlist_lattice.cpp.

В следующем примере мы попытаемся использовать тип<gsl_vector>изGSL(Научная библиотека GNU) в качестве государственного типа в одейте. Мы реализуем это, реализуя обертку вокруг gsl_vector, которая заботится о строительстве / разрушении. Кроме того,Boost.Rangeрасширен таким образом, что он работает с<gsl_vector>s, что также требует реализации нового<gsl_iterator>.

	Note
	odeint уже включает в себя весь код, представленный здесь, см.gsl_wrapper.hpp, поэтому<gsl_vector>s можно использовать прямо из коробки. Следующее описание предназначено только для образовательных целей.

GSL — это библиотека C, поэтому<gsl_vector>не имеет ни конструктора, ни деструктора, ни какой-либо функции<begin>или<end>, ни итераторов вообще. Поэтому, чтобы заставить его работать с odeint, нужно реализовать множество вещей. Обратите внимание, что все работы, показанные здесь, уже включены в odeint, поэтому использование<gsl_vector>s в odeint не требует каких-либо дальнейших корректировок. Мы приводим его в качестве учебного примера. Начнем с определения подходящих конструкторов и деструкторов. Для этого он и создал<state_wrapper><gsl_vector>. Государственные обертки используются степперами для создания и управления временными экземплярами государственных типов:

template<>
struct state_wrapper< gsl_vector* >
{
    typedef double value_type;
    typedef gsl_vector* state_type;
    typedef state_wrapper< gsl_vector* > state_wrapper_type;
    state_type m_v;
    state_wrapper( )
    {
        m_v = gsl_vector_alloc( 1 );
    }
    state_wrapper( const state_wrapper_type &x )
    {
        resize( m_v , x.m_v );
        gsl_vector_memcpy( m_v , x.m_v );
    }
    ~state_wrapper()
    {
        gsl_vector_free( m_v );
    }
};

Эта специализация<state_wrapper>рассказывает о том, как создаются, копируются и уничтожаются gsl_векторы. Далее нам нужно изменить размер, это необходимо, потому что gsl_векторы являются объектами динамических размеров:

template<>
struct is_resizeable< gsl_vector* >
{
    typedef boost::true_type type;
    const static bool value = type::value;
};
template <>
struct same_size_impl< gsl_vector* , gsl_vector* >
{
    static bool same_size( const gsl_vector* x , const gsl_vector* y )
    {
        return x->size == y->size;
    }
};
template <>
struct resize_impl< gsl_vector* , gsl_vector* >
{
    static void resize( gsl_vector* x , const gsl_vector* y )
    {
        gsl_vector_free( x );
        x = gsl_vector_alloc( y->size );
    }
};

До сих пор мы определяли создание / разрушение и изменение размера, но gsl_векторы также не поддерживают итераторы, поэтому мы сначала реализуем итератор gsl:

/*
 * defines an iterator for gsl_vector
 */
class gsl_vector_iterator
      : public boost::iterator_facade< gsl_vector_iterator , double ,
                                       boost::random_access_traversal_tag >
{
public :
    gsl_vector_iterator( void ): m_p(0) , m_stride( 0 ) { }
    explicit gsl_vector_iterator( gsl_vector *p ) : m_p( p->data ) , m_stride( p->stride ) { }
    friend gsl_vector_iterator end_iterator( gsl_vector * );
private :
    friend class boost::iterator_core_access;
    friend class const_gsl_vector_iterator;
    void increment( void ) { m_p += m_stride; }
    void decrement( void ) { m_p -= m_stride; }
    void advance( ptrdiff_t n ) { m_p += n*m_stride; }
    bool equal( const gsl_vector_iterator &other ) const { return this->m_p == other.m_p; }
    bool equal( const const_gsl_vector_iterator &other ) const;
    double& dereference( void ) const { return *m_p; }
    double *m_p;
    size_t m_stride;
};

Аналогичный класс существует для<const>версии итератора. Затем у нас есть функция возврата конечного итератора (аналогично<const>снова):

gsl_vector_iterator end_iterator( gsl_vector *x )
{
    gsl_vector_iterator iter( x );
    iter.m_p += iter.m_stride * x->size;
    return iter;
}

Наконец, добавлены привязки дляBoost.Range:

// template<>
inline gsl_vector_iterator range_begin( gsl_vector *x )
{
    return gsl_vector_iterator( x );
}
// template<>
inline gsl_vector_iterator range_end( gsl_vector *x )
{
    return end_iterator( x );
}

Похожие определения для версий<const>. Это в конечном итоге заставляет одеинт работать с векторами gsl в качестве типов состояний. Полный код для этих связываний находится вgsl_wrapper.hpp. Это может показаться довольно сложным, но имейте в виду, что gsl - это предварительно собранная библиотека C.

Как видно выше, стандартный способ выполнения алгебраических операций на контейнероподобных типах состояний в одейте состоит в том, чтобы итерировать элементы контейнера и выполнять элементы операций по базовому типу значений. Это реализуется с помощью<range_algebra>, который используетBoost.Rangeдля получения итераторов государственных типов. Однако существуют и другие способы осуществления алгебраических операций на контейнерах, одним из которых является определение операторов сложения/умножения для контейнеров непосредственно и затем использование<vector_space_algebra>. При использовании этой алгебры для типа state_type необходимо определить следующие операторы:

Определение этих операторов заставляет ваш тип состояния работать с любым базовым шагомером Рунге-Кутта. Однако, если вы хотите использовать пошаговое управление, требуется больше функциональности. В частности, должны быть выполнены такие операции, какmax_i( |err_i| / (alpha * |s_i|)).errиsявляются типами состояния, альфа является скаляром. Как видите, для получения максимального значения нам нужны элемент мудрого абсолютного значения и деления, а также операция уменьшения. Поэтому для контролируемых шаговых шагов должны быть реализованы следующие вещи:

Здесь мы показываем, как реализовать требуемых операторов на государственном типе. В качестве примера мы определяем новый класс<point3D>, представляющий трехмерный вектор с компонентами x,y,z и определяем для него операторы сложения и скалярного умножения. Мы используем. Операторыуменьшают количество кода, который должен быть написан. Класс для типа точки выглядит следующим образом:

class point3D :
    boost::additive1< point3D ,
    boost::additive2< point3D , double ,
    boost::multiplicative2< point3D , double > > >
{
public:
    double x , y , z;
    point3D()
        : x( 0.0 ) , y( 0.0 ) , z( 0.0 )
    { }
    point3D( const double val )
        : x( val ) , y( val ) , z( val )
    { }
    point3D( const double _x , const double _y , const double _z  )
        : x( _x ) , y( _y ) , z( _z )
    { }
    point3D& operator+=( const point3D &p )
    {
        x += p.x; y += p.y; z += p.z;
        return *this;
    }
    point3D& operator*=( const double a )
    {
        x *= a; y *= a; z *= a;
        return *this;
    }
};

В результатерост. Операторыклассов не должны определять операторов внешнего класса, как<operator+(point3D, point3D)>, потому что об этом заботится библиотека операторов. Обратите внимание, что для простых схем Рунге-Кутта (например,<runge_kutta4>) требуются только операторы<+>и<*>. Если, однако, используется контролируемый степпер, необходимо также указать оператора деления</>, поскольку вычисление термина ошибки включает элемент мудрого деления типов состояния. Кроме того, контролируемые степперы требуют функции<abs>, вычисляющей абсолютное значение элемента для типа состояния:

// only required for steppers with error control
point3D operator/( const point3D &p1 , const point3D &p2 )
{
    return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p2.z );
}
point3D abs( const point3D &p )
{
    return point3D( std::abs(p.x) , std::abs(p.y) , std::abs(p.z) );
}

Наконец, мы должны предоставить специализацию для расчета бесконечности нормы состояния:

// also only for steppers with error control
namespace boost { namespace numeric { namespace odeint {
template<>
struct vector_space_norm_inf< point3D >
{
    typedef double result_type;
    double operator()( const point3D &p ) const
    {
        using std::max;
        using std::abs;
        return max( max( abs( p.x ) , abs( p.y ) ) , abs( p.z ) );
    }
};
} } }

Опять же, обратите внимание, что два последних шага были необходимы только в том случае, если вы хотите использовать контролируемые шаги. Для простых степперов достаточно простого<+=>и<*=>операторов. Определив такой тип точки, мы можем легко выполнить интеграцию в системе Лоренца, явно настраивая<vector_space_algebra>в списке аргументов шаблона степпера:

const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
void lorenz( const point3D &x , point3D &dxdt , const double t )
{
    dxdt.x = sigma * ( x.y - x.x );
    dxdt.y = R * x.x - x.y - x.x * x.z;
    dxdt.z = -b * x.z + x.x * x.y;
}
using namespace boost::numeric::odeint;
int main()
{
    point3D x( 10.0 , 5.0 , 5.0 );
    // point type defines it's own operators -> use vector_space_algebra !
    typedef runge_kutta_dopri5< point3D , double , point3D ,
                                double , vector_space_algebra > stepper;
    int steps = integrate_adaptive( make_controlled<stepper>( 1E-10 , 1E-10 ) , lorenz , x ,
                                    0.0 , 10.0 , 0.1 );
    std::cout << x << std::endl;
    std::cout << "steps: " << steps << std::endl;
}

The whole example can be found in lorenz_point.cpp

	Note
	Для большинства<state_types>одеинт способен автоматически определять правильную алгебру и операции. Но если вы хотите использовать свой<state_type>, как в этом примере с<point3D>, вы должны вручную настроить правильную алгебру / операции, если только ваш<state_type>не работает с выбором по умолчанию<range_algebra>и<default_operations>.