1 /*
2 Copyright 2010-2012 Karsten Ahnert
3 Copyright 2011-2013 Mario Mulansky
4 Copyright 2013 Pascal Germroth
5
6 Distributed under the Boost Software License, Version 1.0.
7 (See accompanying file LICENSE_1_0.txt or
8 copy at http://www.boost.org/LICENSE_1_0.txt)
9 */
10
11
12 #include <iostream>
13 #include <vector>
14
15 #include <boost/numeric/odeint.hpp>
16
17
18
19 //[ rhs_function
20 /* The type of container used to hold the state vector */
21 typedef std::vector< double > state_type;
22
23 const double gam = 0.15;
24
25 /* The rhs of x' = f(x) */
harmonic_oscillator(const state_type & x,state_type & dxdt,const double)26 void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
27 {
28 dxdt[0] = x[1];
29 dxdt[1] = -x[0] - gam*x[1];
30 }
31 //]
32
33
34
35
36
37 //[ rhs_class
38 /* The rhs of x' = f(x) defined as a class */
39 class harm_osc {
40
41 double m_gam;
42
43 public:
harm_osc(double gam)44 harm_osc( double gam ) : m_gam(gam) { }
45
operator ()(const state_type & x,state_type & dxdt,const double)46 void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
47 {
48 dxdt[0] = x[1];
49 dxdt[1] = -x[0] - m_gam*x[1];
50 }
51 };
52 //]
53
54
55
56
57
58 //[ integrate_observer
59 struct push_back_state_and_time
60 {
61 std::vector< state_type >& m_states;
62 std::vector< double >& m_times;
63
push_back_state_and_timepush_back_state_and_time64 push_back_state_and_time( std::vector< state_type > &states , std::vector< double > × )
65 : m_states( states ) , m_times( times ) { }
66
operator ()push_back_state_and_time67 void operator()( const state_type &x , double t )
68 {
69 m_states.push_back( x );
70 m_times.push_back( t );
71 }
72 };
73 //]
74
75 struct write_state
76 {
operator ()write_state77 void operator()( const state_type &x ) const
78 {
79 std::cout << x[0] << "\t" << x[1] << "\n";
80 }
81 };
82
83
main(int,char **)84 int main(int /* argc */ , char** /* argv */ )
85 {
86 using namespace std;
87 using namespace boost::numeric::odeint;
88
89
90 //[ state_initialization
91 state_type x(2);
92 x[0] = 1.0; // start at x=1.0, p=0.0
93 x[1] = 0.0;
94 //]
95
96
97
98 //[ integration
99 size_t steps = integrate( harmonic_oscillator ,
100 x , 0.0 , 10.0 , 0.1 );
101 //]
102
103
104
105 //[ integration_class
106 harm_osc ho(0.15);
107 steps = integrate( ho ,
108 x , 0.0 , 10.0 , 0.1 );
109 //]
110
111
112
113
114
115 //[ integrate_observ
116 vector<state_type> x_vec;
117 vector<double> times;
118
119 steps = integrate( harmonic_oscillator ,
120 x , 0.0 , 10.0 , 0.1 ,
121 push_back_state_and_time( x_vec , times ) );
122
123 /* output */
124 for( size_t i=0; i<=steps; i++ )
125 {
126 cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
127 }
128 //]
129
130
131
132
133
134
135
136 //[ define_const_stepper
137 runge_kutta4< state_type > stepper;
138 integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
139 //]
140
141
142
143
144 //[ integrate_const_loop
145 const double dt = 0.01;
146 for( double t=0.0 ; t<10.0 ; t+= dt )
147 stepper.do_step( harmonic_oscillator , x , t , dt );
148 //]
149
150
151
152
153 //[ define_adapt_stepper
154 typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
155 //]
156
157
158
159 //[ integrate_adapt
160 typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
161 controlled_stepper_type controlled_stepper;
162 integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
163 //]
164
165 {
166 //[integrate_adapt_full
167 double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
168 controlled_stepper_type controlled_stepper(
169 default_error_checker< double , range_algebra , default_operations >( abs_err , rel_err , a_x , a_dxdt ) );
170 integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
171 //]
172 }
173
174
175 //[integrate_adapt_make_controlled
176 integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) ,
177 harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
178 //]
179
180
181
182
183 //[integrate_adapt_make_controlled_alternative
184 integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) ,
185 harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
186 //]
187
188 #ifdef BOOST_NUMERIC_ODEINT_CXX11
189 //[ define_const_stepper_cpp11
190 {
191 runge_kutta4< state_type > stepper;
192 integrate_const( stepper , []( const state_type &x , state_type &dxdt , double t ) {
193 dxdt[0] = x[1]; dxdt[1] = -x[0] - gam*x[1]; }
194 , x , 0.0 , 10.0 , 0.01 );
195 }
196 //]
197
198
199
200 //[ harm_iterator_const_step]
201 std::for_each( make_const_step_time_iterator_begin( stepper , harmonic_oscillator, x , 0.0 , 0.1 , 10.0 ) ,
202 make_const_step_time_iterator_end( stepper , harmonic_oscillator, x ) ,
203 []( std::pair< const state_type & , const double & > x ) {
204 cout << x.second << " " << x.first[0] << " " << x.first[1] << "\n"; } );
205 //]
206 #endif
207
208
209
210
211 }
212