В этом примере мы представляем применение одеина для исследования свойств хаотических детерминистических систем. В математических терминах хаотичный относится к экспоненциальному росту возмущений δ x. Чтобы наблюдать этот экспоненциальный рост, обычно решает уравнения для тангенциальной динамики, которая снова является обычным дифференциальным уравнением. Эти уравнения линейны, но зависят от времени и могут быть получены через

где J - якобиец рассматриваемой системы. δ x также можно интерпретировать как возмущение оригинальной системы. В принципе n этих возмущений существуют, они образуют гиперкуб и развиваются в то время. Экспоненты Lyapunov затем определяются как логарифмические темпы роста возмущений. Если один лиапуновский экспонент больше, то нулевые близлежащие траектории расходятся экспоненциально, следовательно, они хаотичны. Если самый большой экспонент Lyapunov равен нулю, то обычно он сталкивается с периодическим движением. В случае с крупнейшим лиапуновым экспонентом меньшего тогда нулевого сближения до фиксированной точки ожидается. Больше информации об экспонентах Lyapunov и нелинейных динамических системах можно найти во многих учебниках, см., например: E. Ott «Chaos is Dynamical Systems», Cambridge.

Для расчета лиапуновских экспонентов численно одно обычно решает уравнения движения для n пертурбаций и ортонормализует их каждые k шагов. Лиапуновский экспонент является средним логарифмом фактора растяжения каждого возмущения.

Чтобы продемонстрировать, как можно использовать одеинт для определения экспонентов Lyapunov, мы выбираем систему Lorenz. Это одна из самых изученных динамических систем в нелинейном динамическом сообществе. Для стандартных параметров он обладает странным аттрактором с неинтеграционным измерением. Экспоненты Lyapunov принимают значения примерно 0,9, 0 и -12.

Внедрение системы Лоренца

const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
typedef boost::array< double , 3 > lorenz_state_type;
void lorenz( const lorenz_state_type &x , lorenz_state_type &dxdt , double t )
{
    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];
}

Нам также нужно интегрировать набор возмущений. Это делается параллельно оригинальной системе, следовательно, в рамках одной функции системы. Конечно, мы хотим использовать вышеупомянутое определение системы Лоренца, поэтому определение функции системы, включая саму систему Лоренца и возмущение, может выглядеть так:

const size_t n = 3;
const size_t num_of_lyap = 3;
const size_t N = n + n*num_of_lyap;
typedef std::tr1::array< double , N > state_type;
typedef std::tr1::array< double , num_of_lyap > lyap_type;
void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
{
    lorenz( x , dxdt , t );
    for( size_t l=0 ; l<num_of_lyap ; ++l )
    {
        const double *pert = x.begin() + 3 + l * 3;
        double *dpert = dxdt.begin() + 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];
    }
}

Пертурбации хранятся линейно в state_type за состоянием системы Лоренца. Проблема lorenz() и lorenz_with_lyap() с различными типами состояний может быть решена, помещая систему Lorenz внутрь фанктора с искаженными аргументами:

struct lorenz
{
    template< class StateIn , class StateOut , class Value >
    void operator()( const StateIn &x , StateOut &dxdt , Value t )
    {
        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];
    }
};
void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
{
    lorenz()( x , dxdt , t );
    ...
}

Это работает отлично и lorenz_with_lyap может использоваться, например, через

state_type x;
// initialize x..
explicit_rk4< state_type > rk4;
integrate_n_steps( rk4 , lorenz_with_lyap , x , 0.0 , 0.01 , 1000 );

Этот фрагмент кода выполняет 1000 шагов с постоянным размером шага 0,01.

Реальный пример использования в мире для расчета лиапуновых экспонентов системы Лоренца всегда будет включать в себя некоторые переходные шаги, просто для того, чтобы нынешнее состояние лежало на аттракторе, следовательно, оно будет выглядеть как

state_type x;
// initialize x
explicit_rk4< state_type > rk4;
integrate_n_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 );

Проблема в том, что x является полным состоянием, содержащим также возмущения и integrate_n_steps не знает, что он должен использовать только 3 элемента. Подробно олен и его степарды определяют длину рассматриваемой системы путем определения длины состояния. В классических решениях, например, из нумерических рецептов, проблема была решена указателем на состояние и соответствующей длиной, что-то похожее на

void lorenz( double* x , double *dxdt , double t, void* params )
{
    ...
}
int system_length = 3;
rk4( x , system_length , t , dt , lorenz );

Но odeint поддерживает аналогичную и гораздо более сложную концепцию: Boost.Range. Чтобы сделать шагеры и систему готовыми к работе с Boost.Range систему необходимо изменить:

struct lorenz
{
    template< class State , class Deriv >
    void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
    {
        typename boost::range_iterator< const State >::type x = boost::begin( x_ );
        typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
        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];
    }
};

Это в принципе все. Теперь нам нужно только позвонить integrate_n_steps с диапазоном, включая только первые 3 компонента x:

// explicitly choose range_algebra to override default choice of array_algebra
runge_kutta4< state_type , double , state_type , double , range_algebra > rk4;
// perform 10000 transient steps
integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );

	Note
	Обратите внимание, что при использовании Boost.Range, мы должны четко настроить ступенчатый для использования range_algebra, поскольку в противном случае odeint автоматически выбирает array_algebra, что несовместимо с использованием Boost.Range, потому что оригинальный тип состояния является array.

Интегрировав достаточное количество переходных шагов, мы теперь можем вычислить лиапуновых экспонентов:

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 % 100000) )
    {
        cout << t;
        for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
        cout << endl;
    }
}

Полный код можно найти здесь: chaotic_system.cpp