1 /*
2 * Copyright 2011-2012 Mario Mulansky
3 * Copyright 2012-2013 Karsten Ahnert
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 * Example for self defined vector type.
10 */
11
12 #include <vector>
13
14 #include <boost/numeric/odeint.hpp>
15
16 //[my_vector
17 template< size_t MAX_N >
18 class my_vector
19 {
20 typedef std::vector< double > vector;
21
22 public:
23 typedef vector::iterator iterator;
24 typedef vector::const_iterator const_iterator;
25
26 public:
my_vector(const size_t N)27 my_vector( const size_t N )
28 : m_v( N )
29 {
30 m_v.reserve( MAX_N );
31 }
32
my_vector()33 my_vector()
34 : m_v()
35 {
36 m_v.reserve( MAX_N );
37 }
38
39 // ... [ implement container interface ]
40 //]
operator [](const size_t n) const41 const double & operator[]( const size_t n ) const
42 { return m_v[n]; }
43
operator [](const size_t n)44 double & operator[]( const size_t n )
45 { return m_v[n]; }
46
begin()47 iterator begin()
48 { return m_v.begin(); }
49
begin() const50 const_iterator begin() const
51 { return m_v.begin(); }
52
end()53 iterator end()
54 { return m_v.end(); }
55
end() const56 const_iterator end() const
57 { return m_v.end(); }
58
size() const59 size_t size() const
60 { return m_v.size(); }
61
resize(const size_t n)62 void resize( const size_t n )
63 { m_v.resize( n ); }
64
65 private:
66 std::vector< double > m_v;
67
68 };
69
70 //[my_vector_resizeable
71 // define my_vector as resizeable
72
73 namespace boost { namespace numeric { namespace odeint {
74
75 template<size_t N>
76 struct is_resizeable< my_vector<N> >
77 {
78 typedef boost::true_type type;
79 static const bool value = type::value;
80 };
81
82 } } }
83 //]
84
85
86 typedef my_vector<3> state_type;
87
lorenz(const state_type & x,state_type & dxdt,const double t)88 void lorenz( const state_type &x , state_type &dxdt , const double t )
89 {
90 const double sigma( 10.0 );
91 const double R( 28.0 );
92 const double b( 8.0 / 3.0 );
93
94 dxdt[0] = sigma * ( x[1] - x[0] );
95 dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
96 dxdt[2] = -b * x[2] + x[0] * x[1];
97 }
98
99 using namespace boost::numeric::odeint;
100
main()101 int main()
102 {
103 state_type x(3);
104 x[0] = 5.0 ; x[1] = 10.0 ; x[2] = 10.0;
105
106 // make sure resizing is ON
107 BOOST_STATIC_ASSERT( is_resizeable<state_type>::value == true );
108
109 // my_vector works with range_algebra as it implements
110 // the required parts of a container interface
111 // no further work is required
112
113 integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 10.0 , 0.1 );
114 }
115