An important class of ordinary differential equations are so called stiff system which are characterized by two or more time scales of different order. Examples of such systems are found in chemical systems where reaction rates of individual sub-reaction might differ over large ranges, for example:

d S_​1 / dt = - 101 S_​2 - 100 S_​1

d S_​2 / dt = S_​1

In order to efficiently solve stiff systems numerically the Jacobian

J = d f_​i / d x_​j

is needed. Here is the definition of the above example

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

The state type has to be a ublas::vector and the matrix type must by a ublas::matrix since the stiff integrator only accepts these types. However, you might want use non-stiff integrators on this system, too - we will do so later for demonstration. Therefore we want to use the same function also with other state_types, realized by templatizing the 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 )
    {
        ...
    }
};

Now you can use stiff_system in combination with std::vector or boost::array. In the example the explicit time derivative of f(x,t) is introduced separately in the Jacobian. If df / dt = 0 simply fill dfdt with zeros.

A well know solver for stiff systems is the Rosenbrock method. It has a step size control and dense output facilities and can be used like all the other steppers:

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" );

During the integration 71 steps have been done. Comparing to a classical Runge-Kutta solver this is a very good result. For example the Dormand-Prince 5 method with step size control and dense output yields 1531 steps.

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" );

Note, that we have used Boost.Phoenix, a great functional programming library, to create and compose the observer.

The full example can be found here: stiff_system.cpp