Важный класс обыкновенных дифференциальных уравнений называется жесткой системой, которая характеризуется двумя или более временными шкалами разного порядка. Примеры таких систем можно найти в химических системах, где скорость реакции отдельных подреакций может отличаться в больших диапазонах, например:

Для эффективного решения жестких систем численно якобиец

Это необходимо. Вот определение приведенного выше примера

typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
struct stiff_system
{
    void operator()( const vector_type &x , vector_type &dxdt , double /* t */ )
    {
        dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 * x[ 1 ];
        dxdt[ 1 ] = x[ 0 ];
    }
};
struct stiff_system_jacobi
{
    void operator()( const vector_type & /* x */ , matrix_type &J , const double & /* t */ , vector_type &dfdt )
    {
        J( 0 , 0 ) = -101.0;
        J( 0 , 1 ) = -100.0;
        J( 1 , 0 ) = 1.0;
        J( 1 , 1 ) = 0.0;
        dfdt[0] = 0.0;
        dfdt[1] = 0.0;
    }
};

ublas::matrix так как жесткий интегратор принимает только эти типы. Тем не менее, вы можете хотеть использовать интеграторы без тисков на этой системе тоже - мы сделаем это позже для демонстрации. Поэтому мы хотим использовать ту же функцию и с другими типами состояния, реализованными operator():

typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
struct stiff_system
{
    template< class State >
    void operator()( const State &x , State &dxdt , double t )
    {
        ...
    }
};
struct stiff_system_jacobi
{
    template< class State , class Matrix >
    void operator()( const State &x , Matrix &J , const double &t , State &dfdt )
    {
        ...
    }
};

Теперь вы можете использовать<stiff_system>в сочетании с<std::vector>или<boost::array>. В примере явная производная времениf(x,t)вводится отдельно в Якобиане. Еслиdf/dt = 0, просто заполните<dfdt>нулями.

Хорошо известный решатель для жестких систем — метод Розенброка. Он имеет контроль размера шага и плотные производственные мощности и может использоваться, как и все другие степперы:

vector_type x( 2 , 1.0 );
size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< double > >( 1.0e-6 , 1.0e-6 ) ,
        make_pair( stiff_system() , stiff_system_jacobi() ) ,
        x , 0.0 , 50.0 , 0.01 ,
        cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );

За время интеграции было сделано 71 шаг. По сравнению с классическим решателем Рунге-Кутта это очень хороший результат. Например, метод Dormand-Prince 5 с контролем размера шага и плотным выходом дает 1531 шаг.

vector_type x2( 2 , 1.0 );
size_t num_of_steps2 = integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) ,
        stiff_system() , x2 , 0.0 , 50.0 , 0.01 ,
        cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );

Обратите внимание, что мы использовали Boost.Phoenix, большую библиотеку функционального программирования, для создания и составления наблюдателя.