1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
/*
[auto_generated]
libs/numeric/odeint/test/rosenbrock4.cpp
[begin_description]
This file tests the Rosenbrock 4 stepper and its controller and dense output stepper.
[end_description]
Copyright 2011-2012 Karsten Ahnert
Copyright 2011 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)
*/
// disable checked iterator warning for msvc
#include <boost/config.hpp>
#ifdef BOOST_MSVC
#pragma warning(disable:4996)
#endif
#define BOOST_TEST_MODULE odeint_rosenbrock4
#include <utility>
#include <iostream>
#include <boost/test/unit_test.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
#include <boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
typedef double value_type;
typedef boost::numeric::ublas::vector< value_type > state_type;
typedef boost::numeric::ublas::matrix< value_type > matrix_type;
struct sys
{
void operator()( const state_type &x , state_type &dxdt , const value_type &t ) const
{
dxdt( 0 ) = x( 0 ) + 2 * x( 1 );
dxdt( 1 ) = x( 1 );
}
};
struct jacobi
{
void operator()( const state_type &x , matrix_type &jacobi , const value_type &t , state_type &dfdt ) const
{
jacobi( 0 , 0 ) = 1;
jacobi( 0 , 1 ) = 2;
jacobi( 1 , 0 ) = 0;
jacobi( 1 , 1 ) = 1;
dfdt( 0 ) = 0.0;
dfdt( 1 ) = 0.0;
}
};
BOOST_AUTO_TEST_SUITE( rosenbrock4_test )
BOOST_AUTO_TEST_CASE( test_rosenbrock4_stepper )
{
typedef rosenbrock4< value_type > stepper_type;
stepper_type stepper;
typedef stepper_type::state_type state_type;
typedef stepper_type::value_type stepper_value_type;
typedef stepper_type::deriv_type deriv_type;
typedef stepper_type::time_type time_type;
state_type x( 2 ) , xerr( 2 );
x(0) = 0.0; x(1) = 1.0;
stepper.do_step( std::make_pair( sys() , jacobi() ) , x , 0.0 , 0.1 , xerr );
stepper.do_step( std::make_pair( sys() , jacobi() ) , x , 0.0 , 0.1 );
// using std::abs;
// value_type eps = 1E-12;
//
// // compare with analytic solution of above system
// BOOST_CHECK_SMALL( abs( x(0) - 20.0/81.0 ) , eps );
// BOOST_CHECK_SMALL( abs( x(1) - 10.0/9.0 ) , eps );
}
BOOST_AUTO_TEST_CASE( test_rosenbrock4_controller )
{
typedef rosenbrock4_controller< rosenbrock4< value_type > > stepper_type;
stepper_type stepper;
typedef stepper_type::state_type state_type;
typedef stepper_type::value_type stepper_value_type;
typedef stepper_type::deriv_type deriv_type;
typedef stepper_type::time_type time_type;
state_type x( 2 );
x( 0 ) = 0.0 ; x(1) = 1.0;
value_type t = 0.0 , dt = 0.01;
stepper.try_step( std::make_pair( sys() , jacobi() ) , x , t , dt );
}
BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output )
{
typedef rosenbrock4_dense_output< rosenbrock4_controller< rosenbrock4< value_type > > > stepper_type;
typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type;
controlled_stepper_type c_stepper;
stepper_type stepper( c_stepper );
typedef stepper_type::state_type state_type;
typedef stepper_type::value_type stepper_value_type;
typedef stepper_type::deriv_type deriv_type;
typedef stepper_type::time_type time_type;
state_type x( 2 );
x( 0 ) = 0.0 ; x(1) = 1.0;
stepper.initialize( x , 0.0 , 0.1 );
std::pair< value_type , value_type > tr = stepper.do_step( std::make_pair( sys() , jacobi() ) );
stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
}
class rosenbrock4_controller_max_dt_adaptable : public rosenbrock4_controller< rosenbrock4< value_type > >
{
public:
void set_max_dt(value_type max_dt)
{
m_max_dt = max_dt;
}
};
BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output_ref )
{
typedef rosenbrock4_dense_output< boost::reference_wrapper< rosenbrock4_controller_max_dt_adaptable > > stepper_type;
rosenbrock4_controller_max_dt_adaptable c_stepper;
stepper_type stepper( boost::ref( c_stepper ) );
typedef stepper_type::state_type state_type;
typedef stepper_type::value_type stepper_value_type;
typedef stepper_type::deriv_type deriv_type;
typedef stepper_type::time_type time_type;
state_type x( 2 );
x( 0 ) = 0.0 ; x(1) = 1.0;
stepper.initialize( x , 0.0 , 0.1 );
std::pair< value_type , value_type > tr = stepper.do_step( std::make_pair( sys() , jacobi() ) );
stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
// adapt the maximal step size to a small value (smaller than the step size) and make a step
const double max_dt = 1e-8;
c_stepper.set_max_dt(max_dt);
stepper.do_step( std::make_pair( sys() , jacobi() ) );
// check if the step was done with the requested maximal step size
BOOST_CHECK_CLOSE(max_dt, stepper.current_time_step(), 1e-14);
}
BOOST_AUTO_TEST_CASE( test_rosenbrock4_copy_dense_output )
{
typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type;
typedef rosenbrock4_dense_output< controlled_stepper_type > stepper_type;
controlled_stepper_type c_stepper;
stepper_type stepper( c_stepper );
stepper_type stepper2( stepper );
}
BOOST_AUTO_TEST_SUITE_END()
|
