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 <iostream>
16
17 #include <gmpxx.h>
18
19 #include <boost/test/unit_test.hpp>
20 #include <boost/array.hpp>
21
22 #include <boost/mpl/vector.hpp>
23
24 #include <boost/numeric/odeint.hpp>
25 //#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
26
27 using namespace boost::unit_test;
28 using namespace boost::numeric::odeint;
29
30 namespace mpl = boost::mpl;
31
32 const int precision = 1024;
33
34 typedef mpf_class value_type;
35 typedef mpf_class state_type;
36
37 //provide min, max and pow functions for mpf types - required for controlled steppers
min(const value_type a,const value_type b)38 value_type min( const value_type a , const value_type b )
39 {
40 if( a<b ) return a;
41 else return b;
42 }
max(const value_type a,const value_type b)43 value_type max( const value_type a , const value_type b )
44 {
45 if( a>b ) return a;
46 else return b;
47 }
pow(const value_type a,const value_type b)48 value_type pow( const value_type a , const value_type b )
49 {
50 // do calculation in double precision
51 return value_type( std::pow( a.get_d() , b.get_d() ) );
52 }
53
54
55 //provide vector_space reduce:
56
57 namespace boost { namespace numeric { namespace odeint {
58
59 template<>
60 struct vector_space_reduce< state_type >
61 {
62 template< class Op >
operator ()boost::numeric::odeint::vector_space_reduce63 state_type operator()( state_type x , Op op , state_type init ) const
64 {
65 init = op( init , x );
66 return init;
67 }
68 };
69
70 } } }
71
72
constant_system(const state_type & x,state_type & dxdt,value_type t)73 void constant_system( const state_type &x , state_type &dxdt , value_type t )
74 {
75 dxdt = value_type( 1.0 , precision );
76 }
77
78
79 /* check runge kutta stepers */
80 typedef mpl::vector<
81 euler< state_type , value_type , state_type , value_type , vector_space_algebra > ,
82 modified_midpoint< state_type , value_type , state_type , value_type , vector_space_algebra > ,
83 runge_kutta4< state_type , value_type , state_type , value_type , vector_space_algebra > ,
84 runge_kutta4_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
85 runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
86 runge_kutta_cash_karp54< state_type , value_type , state_type , value_type , vector_space_algebra > ,
87 runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > ,
88 runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra >
89 > stepper_types;
90
91
92 template< class Stepper >
93 struct perform_runge_kutta_test {
94
operator ()perform_runge_kutta_test95 void operator()( void )
96 {
97 /* We have to specify the desired precision in advance! */
98 mpf_set_default_prec( precision );
99
100 mpf_t eps_ , unity;
101 mpf_init( eps_ ); mpf_init( unity );
102 mpf_set_d( unity , 1.0 );
103 mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
104 value_type eps( eps_ );
105
106 Stepper stepper;
107 state_type x;
108 x = 0.0;
109
110 stepper.do_step( constant_system , x , 0.0 , 0.1 );
111
112 BOOST_MESSAGE( eps );
113 BOOST_CHECK_MESSAGE( abs( x - value_type( 0.1 , precision ) ) < eps , x - 0.1 );
114 }
115 };
116
117
BOOST_AUTO_TEST_CASE_TEMPLATE(runge_kutta_stepper_test,Stepper,stepper_types)118 BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_stepper_test , Stepper , stepper_types )
119 {
120 perform_runge_kutta_test< Stepper > tester;
121 tester();
122 }
123
124
125 /* check controlled steppers */
126 typedef mpl::vector<
127 controlled_runge_kutta< runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
128 controlled_runge_kutta< runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
129 controlled_runge_kutta< runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
130 bulirsch_stoer< state_type , value_type , state_type , value_type , vector_space_algebra >
131 > controlled_stepper_types;
132
133
134 template< class Stepper >
135 struct perform_controlled_test {
136
operator ()perform_controlled_test137 void operator()( void )
138 {
139 mpf_set_default_prec( precision );
140
141 mpf_t eps_ , unity;
142 mpf_init( eps_ ); mpf_init( unity );
143 mpf_set_d( unity , 1.0 );
144 mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
145 value_type eps( eps_ );
146
147 Stepper stepper;
148 state_type x;
149 x = 0.0;
150
151 value_type t(0.0);
152 value_type dt(0.1);
153
154 stepper.try_step( constant_system , x , t , dt );
155
156 BOOST_MESSAGE( eps );
157 BOOST_CHECK_MESSAGE( abs( x - value_type( 0.1 , precision ) ) < eps , x - 0.1 );
158 }
159 };
160
BOOST_AUTO_TEST_CASE_TEMPLATE(controlled_stepper_test,Stepper,controlled_stepper_types)161 BOOST_AUTO_TEST_CASE_TEMPLATE( controlled_stepper_test , Stepper , controlled_stepper_types )
162 {
163 perform_controlled_test< Stepper > tester;
164 tester();
165 }
166