diff --git a/test/numeric/order_quadrature_formula.cpp b/test/numeric/order_quadrature_formula.cpp index be0b8158..8a87c612 100644 --- a/test/numeric/order_quadrature_formula.cpp +++ b/test/numeric/order_quadrature_formula.cpp @@ -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 - +/* #ifdef BOOST_MSVC #pragma warning(disable:4996) #endif - +*/ #define BOOST_TEST_MODULE order_quadrature_formula #include @@ -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;