odeint implicit euler simple example

Question

I am interested in solving a system of ODEs with odeint library using an implicit scheme and I have difficulties to implement a simple implicit_euler example.

Looking at the documentation, I managed to make work explicit steppers, adaptive ones, as well as the rosenbrock4 stepper. The former seem to be semi-implicit. Therefore I was interested in implementing a fully implicit scheme (and at the same time retrieve the jacobian matrix at each time step). But I did not manage to find documentation and working examples of this stepper. What I have is

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

typedef implicit_euler< double > stepper_IE;
vector_type inout( 2 , 1.0 );
size_t steps = integrate_const( stepper_IE() ,
            std::make_pair( stiff_system() , stiff_system_jacobi() ) ,
            inout , 0.0 , 5.0 , 0.01, streaming_observer( std::cout, x_vec , times ));

The error is the following:

C:\boost_1_55_0\boost\numeric\odeint\stepper\implicit_euler.hpp:94: error: C2064: term does not evaluate to a function taking 3 arguments class does not define an 'operator()' or a user defined conversion operator to a pointer-to-function or reference-to-function that takes appropriate number of arguments

My question now is: does someone know how to make it work, or can someone point me towards a more detailed documentation than this one:

codeproject documentation

or this one:

main odeint page

Thanks

Answer 1

Unfortunately, the implicit Euler method and the Rosenbrock solver (another implicit solver) do not have the same interface. In detail, the implicit Euler expects a function for the Jacobian with this signature

void jacobian( const state_type &x , matrix_type &jacobi , const value_type t );

Hence, you need to change your definition of stiff_system_jacobi to

struct stiff_system_jacobi
{
    void operator()( const vector_type & , matrix_type &J , const double & ) const
    {
        J( 0 , 0 ) = -101.0;
        J( 0 , 1 ) = -100.0;
        J( 1 , 0 ) = 1.0;
        J( 1 , 1 ) = 0.0;
    }
};

If your system is really non autonomous you need to enhance you state type by one additional coordinate which represents the time and has trivial dynamics dt/dt = 1.