boost_mirror/odeint
1
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2025-05-09 23:24:01 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #19 from boostorg/from_headmyshoulder

Browse Source 
Merge from headmyshoulder/odeint-v2
This commit is contained in:
Mario Mulansky 2017-05-09 20:57:34 -07:00 committed by GitHub
commit c0e1cf30f9
15 changed files with 129 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
.travis.yml
View File

@ -19,7 +19,7 @@ addons:
env:
- CXXSTD=""
- CXXSTD="cxxflags='-std=c++0x'"
- CXXSTD="-std=c++0x"
# For now disable gcc on osx as g++4.8 is not yet available
matrix:
@ -28,21 +28,19 @@ matrix:
compiler: gcc
before_install:
# 1.55: http://sourceforge.net/projects/boost/files/boost/1.55.0/boost_1_55_0.tar.bz2/download
# 1.57: http://sourceforge.net/projects/boost/files/boost/1.57.0/boost_1_57_0.tar.bz2/download
# 1.58: http://downloads.sourceforge.net/project/boost/boost/1.58.0/boost_1_58_0.tar.bz2\?r\=\&ts\=1435589970\&use_mirror\=garr
- wget http://downloads.sourceforge.net/project/boost/boost/1.58.0/boost_1_58_0.tar.bz2\?r\=\&ts\=1435589970\&use_mirror\=garr -O /tmp/boost.tar.bz2
# 1.63: https://sourceforge.net/projects/boost/files/boost/1.63.0/boost_1_63_0.tar.bz2/download?use_mirror=superb-dca2
- wget https://sourceforge.net/projects/boost/files/boost/1.63.0/boost_1_63_0.tar.bz2/download?use_mirror=superb-dca2 -O /tmp/boost.tar.bz2
- tar jxf /tmp/boost.tar.bz2
- mv boost_1_58_0 $PWD/boost-trunk
# patch the boost build system - not neccessary with 1.58 anymore
# - patch $PWD/boost-trunk/tools/build/v2/build/toolset.jam toolset.jam.patch
- mv boost_1_63_0 $PWD/boost-trunk
- export BOOST_ROOT="$PWD/boost-trunk"
- cd $BOOST_ROOT
- ./bootstrap.sh
- cd $TRAVIS_BUILD_DIR
- if [ "$TRAVIS_OS_NAME" = "osx" ] && [ "$CC" = "gcc" ]; then export CC=gcc-4.8; fi
- echo $CC
- $CC --version
script:
- $BOOST_ROOT/b2 toolset=$CC $CXXSTD
- $BOOST_ROOT/b2 toolset=$CC cxxflags='$CXXSTD' -j2

2
README.md
View File

@ -1,3 +1,3 @@
[![Bitdeli Badge](https://d2weczhvl823v0.cloudfront.net/headmyshoulder/odeint-v2/trend.png)](https://bitdeli.com/free "Bitdeli Badge")
[![Build Status](https://travis-ci.org/headmyshoulder/odeint-v2.svg?branch=master)](https://travis-ci.org/headmyshoulder/odeint-v2)
odeint is a highly flexible library for solving ordinary differential equations.

1
include/boost/numeric/odeint/integrate/integrate_n_steps.hpp
View File

@ -23,6 +23,7 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/integrate/null_observer.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp>
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
namespace boost {
namespace numeric {

13
include/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
View File

@ -150,6 +150,12 @@ public :
}
void reset(void)
{
m_adams_bashforth.reset();
}
private:
@ -175,7 +181,7 @@ private:
{
m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
m_adams_bashforth.do_step( system , in , t , m_x.m_v , dt );
m_adams_moulton.do_step( system , in , m_x.m_v , t , out , dt , m_adams_bashforth.step_storage() );
m_adams_moulton.do_step( system , in , m_x.m_v , t+dt , out , dt , m_adams_bashforth.step_storage() );
}
else
{
@ -293,6 +299,11 @@ private:
* \param dt The step size.
*/
/**
* \fn adams_bashforth_moulton::reset( void )
* \brief Resets the internal buffers of the stepper.
*/
} // odeint
} // numeric

23
include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
View File

@ -98,6 +98,7 @@ public:
m_interval_sequence( m_k_max+1 ) ,
m_coeff( m_k_max+1 ) ,
m_cost( m_k_max+1 ) ,
m_facmin_table( m_k_max+1 ) ,
m_table( m_k_max ) ,
STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
{
@ -112,21 +113,14 @@ public:
else
m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
m_coeff[i].resize(i);
m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
for( size_t k = 0 ; k < i ; ++k )
{
const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
}
// crude estimate of optimal order
m_current_k_opt = 4;
/* no calculation because log10 might not exist for value_type!
const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
*/
}
reset();
}
@ -282,7 +276,7 @@ public:
{
m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 );
new_h = h_opt[k];
new_h *= m_cost[m_current_k_opt]/m_cost[k];
new_h *= static_cast<value_type>(m_cost[m_current_k_opt])/static_cast<value_type>(m_cost[k]);
} else
new_h = h_opt[m_current_k_opt];
break;
@ -344,6 +338,12 @@ public:
{
m_first = true;
m_last_step_rejected = false;
// crude estimate of optimal order
m_current_k_opt = 4;
/* no calculation because log10 might not exist for value_type!
const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
*/
}
@ -417,7 +417,7 @@ private:
BOOST_USING_STD_MAX();
using std::pow;
value_type expo( 1.0/(2*k+1) );
value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
value_type facmin = m_facmin_table[k];
value_type fac;
if (error == 0.0)
fac=1.0/facmin;
@ -506,6 +506,7 @@ private:
int_vector m_interval_sequence; // stores the successive interval counts
value_matrix m_coeff;
int_vector m_cost; // costs for interval count
value_vector m_facmin_table; // for precomputed facmin to save pow calls
state_table_type m_table; // sequence of states for extrapolation

10
include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
View File

@ -110,6 +110,7 @@ public:
m_interval_sequence( m_k_max+1 ) ,
m_coeff( m_k_max+1 ) ,
m_cost( m_k_max+1 ) ,
m_facmin_table( m_k_max+1 ) ,
m_table( m_k_max ) ,
m_mp_states( m_k_max+1 ) ,
m_derivs( m_k_max+1 ) ,
@ -125,9 +126,13 @@ public:
m_interval_sequence[i] = 2 + 4*i; // 2 6 10 14 ...
m_derivs[i].resize( m_interval_sequence[i] );
if( i == 0 )
{
m_cost[i] = m_interval_sequence[i];
else
} else
{
m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
}
m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) );
m_coeff[i].resize(i);
for( size_t k = 0 ; k < i ; ++k )
{
@ -429,7 +434,7 @@ private:
using std::pow;
value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]);
value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
value_type facmin = m_facmin_table[k];
value_type fac;
if (error == 0.0)
fac = static_cast<value_type>(1)/facmin;
@ -692,6 +697,7 @@ private:
int_vector m_interval_sequence; // stores the successive interval counts
value_matrix m_coeff;
int_vector m_cost; // costs for interval count
value_vector m_facmin_table; // for precomputed facmin to save pow calls
state_vector_type m_table; // sequence of states for extrapolation

2
include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
View File

@ -143,7 +143,7 @@ public:
error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
error);
time_type dt_old = dt;
// time_type dt_old = dt; unused variable warning
//error too small - increase dt and keep the evolution and limit scaling factor to 5.0
dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
pow(error, static_cast<value_type>(-1) / stepper_order);

2
include/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
View File

@ -163,6 +163,8 @@ public:
if( m_max_dt != static_cast<time_type>(0) )
{
dt = detail::min_abs(m_max_dt, dt_new);
} else {
dt = dt_new;
}
m_last_rejected = false;
return success;

7
test/adams_bashforth.cpp
View File

@ -208,8 +208,13 @@ BOOST_AUTO_TEST_CASE( test_auto_initialization )
adams.do_step( lorenz() , x , 0.0 , x , 0.1 );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 2 ) );
adams.reset();
adams.do_step( lorenz() , x , 0.0 , x , 0.1 );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 2 ) );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 3 ) );
adams.do_step( lorenz() , x , 0.0 , x , 0.1 );
BOOST_CHECK_EQUAL( adams.initializing_stepper().do_count , size_t( 4 ) );
}
BOOST_AUTO_TEST_CASE( test_manual_initialization )

4
test/adams_bashforth_moulton.cpp
View File

@ -104,8 +104,8 @@ BOOST_AUTO_TEST_CASE( test_instantiation )
s4.do_step( lorenz() , x , t , dt );
s5.do_step( lorenz() , x , t , dt );
s6.do_step( lorenz() , x , t , dt );
// s7.do_step( lorenz() , x , t , dt );
// s8.do_step( lorenz() , x , t , dt );
// s7.do_step( lorenz() , x , t , dt );
// s8.do_step( lorenz() , x , t , dt );
}

2
test/numeric/order_quadrature_formula.cpp
View File

@ -20,6 +20,8 @@
#include <iostream>
#include <cmath>
#include "boost/format.hpp"
#include <boost/test/unit_test.hpp>
#include <boost/mpl/vector.hpp>

1
test/regression/Jamfile.v2
View File

@ -26,5 +26,6 @@ test-suite "odeint"
[ run regression_147.cpp ]
[ compile regression_149.cpp ]
[ run regression_168.cpp ]
[ run regression_189.cpp ]
: <testing.launcher>valgrind
;

72
test/regression/regression_189.cpp Normal file
View File

@ -0,0 +1,72 @@
/*
[begin_description]
Test case for issue 189:
Controlled Rosenbrock stepper fails to increase step size
[end_description]
Copyright 2016 Karsten Ahnert
Copyright 2016 Mario Mulansky
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)
*/
#define BOOST_TEST_MODULE odeint_regression_189
#include <boost/numeric/odeint.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/phoenix/core.hpp>
#include <boost/phoenix/operator.hpp>
using namespace boost::numeric::odeint;
namespace phoenix = boost::phoenix;
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;
}
};
BOOST_AUTO_TEST_CASE( regression_189 )
{
vector_type x( 2 , 1.0 );
size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< double > >( 1.0e-6 , 1.0e-6 ) ,
std::make_pair( stiff_system() , stiff_system_jacobi() ) ,
x , 0.0 , 50.0 , 0.01 ,
std::cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
// regression: number of steps should be 74
size_t num_of_steps_expected = 74;
BOOST_CHECK_EQUAL( num_of_steps , num_of_steps_expected );
vector_type x2( 2 , 1.0 );
size_t num_of_steps2 = integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) ,
stiff_system() , x2 , 0.0 , 50.0 , 0.01 ,
std::cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
num_of_steps_expected = 1531;
BOOST_CHECK_EQUAL( num_of_steps2 , num_of_steps_expected );
}

1
test/runge_kutta_concepts.cpp
View File

@ -88,6 +88,7 @@ struct perform_stepper_test
typedef T vector_space_type;
void operator()( void ) const
{
using std::abs;
vector_space_type x;
x = 2.0;
Stepper stepper;

2
test/runge_kutta_controlled_concepts.cpp
View File

@ -87,6 +87,7 @@ struct perform_controlled_stepper_test
typedef T vector_space_type;
void operator()( void ) const
{
using std::abs;
vector_space_type x;
x = 2.0;
ControlledStepper controlled_stepper;
@ -134,6 +135,7 @@ struct perform_controlled_stepper_test< ControlledStepper , vector_space_type >
{
void operator()( void ) const
{
using std::abs;
vector_space_type x;
x = 2.0;
ControlledStepper controlled_stepper;