1 | /*
|
---|
2 | Copyright 2009-2012 Karsten Ahnert
|
---|
3 | Copyright 2009-2012 Mario Mulansky
|
---|
4 |
|
---|
5 | Distributed under the Boost Software License, Version 1.0.
|
---|
6 | (See accompanying file LICENSE_1_0.txt or
|
---|
7 | copy at http://www.boost.org/LICENSE_1_0.txt)
|
---|
8 | */
|
---|
9 |
|
---|
10 |
|
---|
11 | #include <iostream>
|
---|
12 | #include <vector>
|
---|
13 |
|
---|
14 | #include <boost/numeric/odeint.hpp>
|
---|
15 |
|
---|
16 |
|
---|
17 |
|
---|
18 | //[ rhs_function
|
---|
19 | /* The type of container used to hold the state vector */
|
---|
20 | typedef std::vector< double > state_type;
|
---|
21 |
|
---|
22 | const double gam = 0.15;
|
---|
23 |
|
---|
24 | /* The rhs of x' = f(x) */
|
---|
25 | void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
|
---|
26 | {
|
---|
27 | dxdt[0] = x[1];
|
---|
28 | dxdt[1] = -x[0] - gam*x[1];
|
---|
29 | }
|
---|
30 |
|
---|
31 |
|
---|
32 | //[ integrate_observer
|
---|
33 | struct push_back_state_and_time
|
---|
34 | {
|
---|
35 | std::vector< state_type >& m_states;
|
---|
36 | std::vector< double >& m_times;
|
---|
37 |
|
---|
38 | push_back_state_and_time( std::vector< state_type > &states , std::vector< double > × )
|
---|
39 | : m_states( states ) , m_times( times ) { }
|
---|
40 |
|
---|
41 | void operator()( const state_type &x , double t )
|
---|
42 | {
|
---|
43 | m_states.push_back( x );
|
---|
44 | m_times.push_back( t );
|
---|
45 | }
|
---|
46 | };
|
---|
47 |
|
---|
48 |
|
---|
49 | int main(int /* argc */ , char** /* argv */ )
|
---|
50 | {
|
---|
51 | using namespace std;
|
---|
52 | using namespace boost::numeric::odeint;
|
---|
53 |
|
---|
54 |
|
---|
55 | //[ state_initialization
|
---|
56 | state_type x(2);
|
---|
57 | x[0] = 1.0; // start at x=1.0, p=0.0
|
---|
58 | x[1] = 0.0;
|
---|
59 |
|
---|
60 |
|
---|
61 | //[ integrate_observ
|
---|
62 | vector<state_type> x_vec;
|
---|
63 | vector<double> times;
|
---|
64 |
|
---|
65 |
|
---|
66 | //[ define_const_stepper
|
---|
67 | bulirsch_stoer_dense_out< state_type > stepper;
|
---|
68 | stepper.adjust_size(x_vec);
|
---|
69 | size_t steps = integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01, push_back_state_and_time( x_vec , times ) );
|
---|
70 | //]
|
---|
71 |
|
---|
72 | for( size_t i=0; i<=steps; i++ )
|
---|
73 | {
|
---|
74 | cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
|
---|
75 | }
|
---|
76 | }
|
---|