odeint/test/numeric/order_quadrature_formula.cpp

/* Boost numeric test for ordersof quadrature formulas test file

 Copyright 2012 Mario Mulansky
 Copyright 2012 Karsten Ahnert

 Distributed under the Boost Software License, Version 1.0.
 (See accompanying file LICENSE_1_0.txt or
 copy at http://www.boost.org/LICENSE_1_0.txt)
*/

// disable checked iterator warning for msvc
#include <boost/config.hpp>

#ifdef BOOST_MSVC
    #pragma warning(disable:4996)
#endif

#define BOOST_TEST_MODULE order_quadrature_formula

#include <iostream>
#include <cmath>

#include <boost/test/unit_test.hpp>

#include <boost/mpl/vector.hpp>

#include <boost/numeric/odeint.hpp>

#include <boost/numeric/ublas/vector.hpp>

using namespace boost::unit_test;
using namespace boost::numeric::odeint;
namespace mpl = boost::mpl;

typedef double value_type;
typedef value_type time_type;
typedef boost::numeric::ublas::vector< double > state_type;

BOOST_AUTO_TEST_SUITE( order_of_convergence_test )

int p;
value_type tolerance = 1.0e-13;

/* generic test for all steppers that support integrate_const */
template< class Stepper >
struct integrate_const_test
{
    double error;
    double maxError;
    void operator()( int nSteps = 1 )
    {
	Stepper stepper;
	std::cout << boost::format( "%-20i%-20i%-30E%-20E\n" )
	    % estimatedOrder( nSteps ) %  definedOrder() % error % maxError;

	BOOST_REQUIRE_EQUAL( estimatedOrder( nSteps ), definedOrder() );
    }

    const int definedOrder(){
	const Stepper stepper;
	return stepper.order();
    }
    /*
    the order of the stepper is estimated by trying to solve the ODE
    x'(t) = t^p
    until the errors are too big to be justified by finite precision.
    the first value p for which the problem is *not* solved with
    precision `tolerance` is the estimate for the order of the scheme.
     */
    const int estimatedOrder( int nSteps = 1 )
    {
	state_type x(1);
	double t;
	maxError = 0;
	for ( p = 0; p < 20; p++ )
	    {
		Stepper stepper;
		x[0] = 0.0;
		t = 0;
		for ( int i = 0; i < nSteps; i++ ){
		    stepper.do_step( rhs, x, t, 1.0/nSteps );
		    t += 1.0/nSteps;
		    error = fabs( x[0] - pow( t, ( 1.0 + p ) )/( 1.0 + p ) );
		    if ( error > tolerance )
			return p;
		    else if ( error > maxError )
			maxError = error;
		}
	    }
	return p;
    }


private:
    static void rhs( const state_type &x , state_type &dxdt , const time_type t )
    {
	dxdt[0] = pow( t, p );
    }
};


template< class Stepper >
struct integrate_const_test_initialize
{
    value_type error;
    value_type maxError;

    void operator()( int nSteps = 16 )
    {
	std::cout << boost::format( "%-20i%-20i%-30E%-20E\n" )
	    % estimatedOrder( nSteps ) %  definedOrder() % error % maxError;

	BOOST_REQUIRE_EQUAL( estimatedOrder( nSteps ), definedOrder() );
    }

    const int definedOrder(){
	const Stepper stepper;
	return stepper.order();
    }
    /*
    just like the other version of estimatedOrder(), but with
    initization performed by the fehlberg stepper ( order 8 )

    if the default initialization ( runge kutta 4 ) is used,
    estimatedOrder will never return an order greater than 4.
    */
    const int estimatedOrder( int nSteps = 16 )
    {
	state_type x(1);
	value_type t;
	maxError = -1.0;
	for ( p = 0; p < 10; p++ )
	    {
		Stepper stepper;
		x[0] = 0.0;
		t = 0;
		// use a high order method to initialize.
		stepper.initialize( runge_kutta_fehlberg78< state_type >(),
					 rhs, x, t, 1.0/nSteps );
		while ( t < 1 ){
		    stepper.do_step( rhs, x, t, 1.0/nSteps );
		    t += 1.0/nSteps;
		    error = fabs( x[0]-pow( t, 1.0 + p)/(1.0+p) );
		    if ( error > tolerance )
			return p;
		    else if ( error > maxError )
			maxError = error;
		}
	    }
	return p;
    }


private:
    static void rhs( const state_type &x , state_type &dxdt , const time_type t )
    {
	dxdt[0] = pow( t, p );
    }
};


typedef mpl::vector<
    euler< state_type > ,
    modified_midpoint< state_type > ,
    runge_kutta4< state_type > ,
    runge_kutta4_classic< state_type > ,
    runge_kutta_cash_karp54_classic< state_type > ,
    runge_kutta_cash_karp54< state_type > ,
    runge_kutta_dopri5< state_type > ,
    runge_kutta_fehlberg78< state_type >
    > runge_kutta_steppers;

typedef mpl::vector<
    adams_bashforth< 2, state_type > ,
    adams_bashforth< 3, state_type > ,
    adams_bashforth< 4, state_type > ,
    adams_bashforth< 5, state_type > ,
    adams_bashforth< 6, state_type > ,
    adams_bashforth< 7, state_type > ,
    adams_bashforth< 8, state_type > ,
    adams_bashforth_moulton< 2, state_type > ,
    adams_bashforth_moulton< 3, state_type > ,
    adams_bashforth_moulton< 4, state_type > ,
    adams_bashforth_moulton< 5, state_type > ,
    adams_bashforth_moulton< 6, state_type > ,
    adams_bashforth_moulton< 7, state_type > ,
    adams_bashforth_moulton< 8, state_type >
    > adams_steppers;

BOOST_AUTO_TEST_CASE( print_header )
{
    std::cout << boost::format( "%-20s%-20s%-30s%-20s\n" )
	% "Estimated order" % "defined order"
	% "first significant error" % "biggest insignificant error";
}

BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers )
{
    integrate_const_test< Stepper > tester;
    tester();
}

BOOST_AUTO_TEST_CASE_TEMPLATE( adams_moultion_test , Stepper, adams_steppers )
{
    integrate_const_test_initialize< Stepper > tester;
    tester();
}

BOOST_AUTO_TEST_SUITE_END()