196 lines
6.8 KiB
C++
196 lines
6.8 KiB
C++
/* Boost numeric test for orders of quadrature formulas test file
|
|
|
|
Copyright 2015 Gregor de Cillia
|
|
Copyright 2015 Mario Mulansky <mario.mulansky@gmx.net>
|
|
|
|
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/format.hpp"
|
|
|
|
#include <boost/test/unit_test.hpp>
|
|
|
|
#include <boost/mpl/vector.hpp>
|
|
|
|
#include <boost/numeric/odeint.hpp>
|
|
|
|
#include <boost/numeric/ublas/vector.hpp>
|
|
|
|
#include <boost/format.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 value_type state_type;
|
|
|
|
BOOST_AUTO_TEST_SUITE( order_of_convergence_test )
|
|
|
|
/* defines the simple monomial f(t) = (p+1) * t^p.*/
|
|
struct monomial
|
|
{
|
|
int power;
|
|
|
|
monomial(int p = 0) : power( p ){};
|
|
|
|
void operator()( const state_type &x , state_type &dxdt , const time_type t )
|
|
{
|
|
dxdt = ( 1.0 + power ) * pow( t, power );
|
|
}
|
|
};
|
|
|
|
|
|
/* generic test for all steppers that support integrate_const */
|
|
template< class Stepper >
|
|
struct stepper_order_test
|
|
{
|
|
void operator()( int steps = 1 )
|
|
{
|
|
const int estimated_order = estimate_order( steps );
|
|
const int defined_order = Stepper::order_value;
|
|
|
|
std::cout << boost::format( "%-20i%-20i\n" )
|
|
% estimated_order % defined_order;
|
|
|
|
BOOST_REQUIRE_EQUAL( estimated_order, defined_order );
|
|
}
|
|
|
|
/*
|
|
the order of the stepper is estimated by trying to solve the ODE
|
|
x'(t) = (p+1) * 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 within the
|
|
finite precision tolerance is the estimate for the order of the scheme.
|
|
*/
|
|
int estimate_order( int steps )
|
|
{
|
|
const double dt = 1.0/steps;
|
|
const double tolerance = steps*1E-15;
|
|
int p;
|
|
for( p = 0; true; p++ )
|
|
{
|
|
// begin with x'(t) = t^0 = 1
|
|
// => x (t) = t
|
|
// then use x'(t) = 2*t^1
|
|
// => x (t) = t^2
|
|
// ...
|
|
state_type x = 0.0;
|
|
|
|
double t = integrate_n_steps( Stepper(), monomial( p ), x, 0.0, dt,
|
|
steps );
|
|
if( fabs( x - pow( t, ( 1.0 + p ) ) ) > tolerance )
|
|
break;
|
|
}
|
|
// the smallest power p for which the test failed is the estimated order,
|
|
// as the solution for this power is x(t) = t^{p+1}
|
|
return 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, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer, runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth< 3, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer, runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth< 4, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer, runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth< 5, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer, runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth< 6, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer, runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth< 7, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer, runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth< 8, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer, runge_kutta_fehlberg78< state_type > >
|
|
> ab_steppers;
|
|
|
|
|
|
typedef mpl::vector<
|
|
adams_bashforth_moulton< 2, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer,
|
|
runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth_moulton< 3, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer,
|
|
runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth_moulton< 4, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer,
|
|
runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth_moulton< 5, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer,
|
|
runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth_moulton< 6, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer,
|
|
runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth_moulton< 7, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer,
|
|
runge_kutta_fehlberg78< state_type > >,
|
|
adams_bashforth_moulton< 8, state_type, double, state_type, double,
|
|
vector_space_algebra, default_operations,
|
|
initially_resizer,
|
|
runge_kutta_fehlberg78< state_type > >
|
|
> abm_steppers;
|
|
|
|
|
|
BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers )
|
|
{
|
|
stepper_order_test< Stepper > tester;
|
|
tester(10);
|
|
}
|
|
|
|
|
|
BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, ab_steppers )
|
|
{
|
|
stepper_order_test< Stepper > tester;
|
|
tester(16);
|
|
}
|
|
|
|
|
|
BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moultion_test , Stepper, abm_steppers )
|
|
{
|
|
stepper_order_test< Stepper > tester;
|
|
tester(16);
|
|
}
|
|
|
|
BOOST_AUTO_TEST_SUITE_END()
|