1 /* Boost check_gmp.cpp test file
2
3 Copyright 2010-2012 Mario Mulansky
4 Copyright 2011-2012 Karsten Ahnert
5
6 This file tests the odeint library with the gmp arbitrary precision types
7
8 Distributed under the Boost Software License, Version 1.0.
9 (See accompanying file LICENSE_1_0.txt or
10 copy at http://www.boost.org/LICENSE_1_0.txt)
11 */
12
13 #define BOOST_TEST_MODULE odeint_gmp
14
15 #include <gmpxx.h>
16
17 #include <boost/test/unit_test.hpp>
18 #include <boost/array.hpp>
19
20 #include <boost/mpl/vector.hpp>
21
22 #include <boost/numeric/odeint.hpp>
23 //#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
24
25 using namespace boost::unit_test;
26 using namespace boost::numeric::odeint;
27
28 namespace mpl = boost::mpl;
29
30 const int precision = 1024;
31
32 typedef mpf_class value_type;
33 typedef mpf_class state_type;
34
35 //provide min, max and pow functions for mpf types - required for controlled steppers
min(const value_type a,const value_type b)36 value_type min( const value_type a , const value_type b )
37 {
38 if( a<b ) return a;
39 else return b;
40 }
max(const value_type a,const value_type b)41 value_type max( const value_type a , const value_type b )
42 {
43 if( a>b ) return a;
44 else return b;
45 }
pow(const value_type a,const value_type b)46 value_type pow( const value_type a , const value_type b )
47 {
48 // do the calculation in double precision...
49 return value_type( std::pow( a.get_d() , b.get_d() ) );
50 }
51
52
53 //provide vector_space reduce:
54
55 namespace boost { namespace numeric { namespace odeint {
56
57 template<>
58 struct vector_space_reduce< state_type >
59 {
60 template< class Op >
operator ()boost::numeric::odeint::vector_space_reduce61 state_type operator()( state_type x , Op op , state_type init ) const
62 {
63 init = op( init , x );
64 return init;
65 }
66 };
67
68 } } }
69
70
constant_system(const state_type & x,state_type & dxdt,value_type t)71 void constant_system( const state_type &x , state_type &dxdt , value_type t )
72 {
73 dxdt = value_type( 1.0 , precision );
74 }
75
76 /* check runge kutta stepers */
77 typedef mpl::vector<
78 euler< state_type , value_type , state_type , value_type , vector_space_algebra > ,
79 modified_midpoint< state_type , value_type , state_type , value_type , vector_space_algebra > ,
80 runge_kutta4< state_type , value_type , state_type , value_type , vector_space_algebra > ,
81 runge_kutta4_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
82 runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
83 runge_kutta_cash_karp54< state_type , value_type , state_type , value_type , vector_space_algebra > ,
84 runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > ,
85 runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra >
86 > stepper_types;
87
88
89 template< class Stepper >
90 struct perform_integrate_const_test {
91
operator ()perform_integrate_const_test92 void operator()( void )
93 {
94 /* We have to specify the desired precision in advance! */
95 mpf_set_default_prec( precision );
96
97 mpf_t eps_ , unity;
98 mpf_init( eps_ ); mpf_init( unity );
99 mpf_set_d( unity , 1.0 );
100 mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
101 value_type eps( eps_ );
102
103 Stepper stepper;
104 state_type x;
105 x = 0.0;
106 value_type t0( 0.0 );
107 value_type tend( 1.0 );
108 value_type dt(0.1);
109
110 integrate_const( stepper , constant_system , x , t0 , tend , dt );
111
112 x = 0.0;
113 t0 = 0.0;
114 dt = 0.1;
115 size_t steps = 10;
116
117 integrate_n_steps( stepper , constant_system , x , t0 , dt , steps );
118
119 BOOST_CHECK_MESSAGE( abs( x - 10*dt ) < eps , x );
120 }
121 };
122
123
BOOST_AUTO_TEST_CASE_TEMPLATE(integrate_const_test,Stepper,stepper_types)124 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_const_test , Stepper , stepper_types )
125 {
126 perform_integrate_const_test< Stepper > tester;
127 tester();
128 }
129
130
131 typedef mpl::vector<
132 controlled_runge_kutta< runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
133 controlled_runge_kutta< runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
134 controlled_runge_kutta< runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
135 bulirsch_stoer< state_type , value_type , state_type , value_type , vector_space_algebra >
136 > controlled_stepper_types;
137
138
139 template< class Stepper >
140 struct perform_integrate_adaptive_test {
141
operator ()perform_integrate_adaptive_test142 void operator()( void )
143 {
144 mpf_set_default_prec( precision );
145
146 mpf_t eps_ , unity;
147 mpf_init( eps_ ); mpf_init( unity );
148 mpf_set_d( unity , 1.0 );
149 mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
150 value_type eps( eps_ );
151
152 Stepper stepper;
153 state_type x;
154 x = 0.0;
155 value_type t0( 0.0 );
156 value_type tend( 1.0 );
157 value_type dt(0.1);
158
159 integrate_adaptive( stepper , constant_system , x , t0 , tend , dt );
160
161 BOOST_CHECK_MESSAGE( abs( x - tend ) < eps , x - 0.1 );
162 }
163 };
164
BOOST_AUTO_TEST_CASE_TEMPLATE(integrate_adaptive__test,Stepper,controlled_stepper_types)165 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_adaptive__test , Stepper , controlled_stepper_types )
166 {
167 perform_integrate_adaptive_test< Stepper > tester;
168 tester();
169 }
170