boost/odeint
2
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2026-02-09 11:22:10 +00:00
Code Issues Projects Releases Wiki Activity

applied recommendations mentioned in #146

Browse Source 
* state_type is now a double
* rhs is global
* p is a member of rhs, called exponent
* 'main loop' is of the form
	do{ exponent++; ... } while ( error < tolerance );
This commit is contained in:
Gregor de Cillia
2015-01-20 01:43:52 +01:00
parent b9fa3d5a1c
commit c7d40ae0f1
1 changed files with 85 additions and 71 deletions
Download Patch File Download Diff File Expand all files Collapse all files

156
test/numeric/order_quadrature_formula.cpp
View File

@@ -1,4 +1,4 @@
/* Boost numeric test for ordersof quadrature formulas test file
/* Boost numeric test for orders of quadrature formulas test file
Copyright 2012 Mario Mulansky
Copyright 2012 Karsten Ahnert
@@ -10,11 +10,11 @@
// 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>
@@ -34,26 +34,36 @@ namespace mpl = boost::mpl;
typedef double value_type;
typedef value_type time_type;
typedef boost::numeric::ublas::vector< double > state_type;
typedef value_type state_type;
BOOST_AUTO_TEST_SUITE( order_of_convergence_test )
int p;
value_type tolerance = 1.0e-13;
struct monome
{
int exponent;
monome() : exponent( 0 ){};
void operator()( const state_type &x , state_type &dxdt , const time_type t ){
dxdt = ( 1.0 + exponent )*pow( t, exponent );
}
};
monome rhs;
/* generic test for all steppers that support integrate_const */
template< class Stepper >
struct integrate_const_test
{
double error;
double maxError;
int estimatedOrder;
void operator()( int nSteps = 1 )
{
Stepper stepper;
std::cout << boost::format( "%-20i%-20i%-30E%-20E\n" )
% estimatedOrder( nSteps ) % definedOrder() % error % maxError;
estimateOrder();
std::cout << boost::format( "%-20i%-20i%-20E\n" )
% estimatedOrder % definedOrder() % maxError;
BOOST_REQUIRE_EQUAL( estimatedOrder( nSteps ), definedOrder() );
BOOST_REQUIRE_EQUAL( estimatedOrder, definedOrder() );
}
const int definedOrder(){
@@ -67,34 +77,33 @@ struct integrate_const_test
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 )
void estimateOrder( int nSteps = 1 )
{
state_type x(1);
state_type x;
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;
}
rhs.exponent = -1;
do{
// begin with x'(t) = ( t^0 )/1 = 1
// => x (t) = t
// then use x'(t) = ( t^1 )/2 = t/2
// => x (t) = t^2
// ...
rhs.exponent++;
Stepper stepper;
x = 0.0;
t = 0;
for ( int i = 0; i < nSteps; i++ ){
stepper.do_step( rhs, x, t, 1.0/nSteps );
t += 1.0/nSteps;
// compute the error using the exact solution x(t)=t^(1+p)
value_type error = fabs( x - pow( t, ( 1.0 + rhs.exponent ) ) );
maxError = fmax( error, maxError );
}
return p;
}
private:
static void rhs( const state_type &x , state_type &dxdt , const time_type t )
{
dxdt[0] = pow( t, p );
}
while ( maxError < tolerance );
// return the first exponent for which the test failed
estimatedOrder = rhs.exponent;
}
};
@@ -102,15 +111,16 @@ private:
template< class Stepper >
struct integrate_const_test_initialize
{
value_type error;
int estimatedOrder;
value_type maxError;
void operator()( int nSteps = 16 )
{
std::cout << boost::format( "%-20i%-20i%-30E%-20E\n" )
% estimatedOrder( nSteps ) % definedOrder() % error % maxError;
estimateOrder( nSteps );
std::cout << boost::format( "%-20i%-20i%-30E\n" )
% estimatedOrder % definedOrder() % maxError;
BOOST_REQUIRE_EQUAL( estimatedOrder( nSteps ), definedOrder() );
BOOST_REQUIRE_EQUAL( estimatedOrder, definedOrder() );
}
const int definedOrder(){
@@ -118,43 +128,39 @@ struct integrate_const_test_initialize
return stepper.order();
}
/*
just like the other version of estimatedOrder(), but with
just like the other version of estimateOrder(), 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.
the estimated order will never be greater than 4.
*/
const int estimatedOrder( int nSteps = 16 )
void estimateOrder( 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;
}
state_type x;
time_type t;
const time_type dt = 1.0/nSteps;
maxError = 0.0;
rhs.exponent = -1;
do{
rhs.exponent++;
// construct the stepper inside the for loop to reset the
// step storage
Stepper stepper;
x = 0.0;
t = 0.0;
// use a high order method to initialize.
stepper.initialize( runge_kutta_fehlberg78< state_type >(),
rhs, x, t, dt );
while ( t < 1 ){
stepper.do_step( rhs, x, t, 1.0/nSteps );
t += 1.0/nSteps;
// compute the error using the exact solution x(t)=t^(1+p)
value_type error = fabs( x - pow( t, 1.0 + rhs.exponent) );
maxError = fmax( error, maxError );
}
return p;
}
private:
static void rhs( const state_type &x , state_type &dxdt , const time_type t )
{
dxdt[0] = pow( t, p );
} while( maxError < tolerance );
// return the first exponent for which the test failed
estimatedOrder = rhs.exponent;
}
};
@@ -189,9 +195,9 @@ typedef mpl::vector<
BOOST_AUTO_TEST_CASE( print_header )
{
std::cout << boost::format( "%-20s%-20s%-30s%-20s\n" )
std::cout << boost::format( "%-20s%-20s%-30s\n" )
% "Estimated order" % "defined order"
% "first significant error" % "biggest insignificant error";
% "maxError";
}
BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers )
@@ -200,6 +206,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers
tester();
}
BOOST_AUTO_TEST_CASE( print_seperator )
{
std::cout << "-------------------"
<< " ABM STEPPERS "
<< "-------------------"
<< std::endl;
}
BOOST_AUTO_TEST_CASE_TEMPLATE( adams_moultion_test , Stepper, adams_steppers )
{
integrate_const_test_initialize< Stepper > tester;