1 [/============================================================================
4 Copyright 2011-2012 Karsten Ahnert
5 Copyright 2011-2013 Mario Mulansky
6 Copyright 2012 Sylwester Arabas
8 Use, modification and distribution is subject to the Boost Software License,
9 Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
10 http://www.boost.org/LICENSE_1_0.txt)
11 =============================================================================/]
13 [section State types, algebras and operations]
15 In odeint the stepper algorithms are implemented independently of the
16 underlying fundamental mathematical operations.
17 This is realized by giving the user full control over the state type and the
18 mathematical operations for this state type.
19 Technically, this is done by introducing three concepts: StateType, Algebra,
21 Most of the steppers in odeint expect three class types fulfilling these
22 concepts as template parameters.
23 Note that these concepts are not fully independent of each other but rather a
24 valid combination must be provided in order to make the steppers work.
25 In the following we will give some examples on reasonable
26 state_type-algebra-operations combinations.
27 For the most common state types, like `vector<double>` or `array<double,N>`
28 the default values range_algebra and default_operations are perfectly fine and
29 odeint can be used as is without worrying about algebra/operations at all.
31 [important state_type, algebra and operations are not independent, a valid
32 combination must be provided to make odeint work properly]
34 Moreover, as odeint handles the memory required for intermediate temporary
35 objects itself, it also needs knowledge about how to create state_type objects
36 and maybe how to allocate memory (resizing).
37 All in all, the following things have to be taken care of when odeint is used
38 with non-standard state types:
40 * construction/destruction
41 * resizing (if possible/required)
42 * algebraic operations
44 Again, odeint already provides basic interfaces for most of the usual state
46 So if you use a `std::vector`, or a `boost::array` as state type no additional
47 work is required, they just work out of the box.
49 [section Construction/Resizing]
51 We distinguish between two basic state types: fixed sized and dynamically
53 For fixed size state types the default constructor `state_type()` already
54 allocates the required memory, prominent example is `boost::array<T,N>`.
55 Dynamically sized types have to be resized to make sure enough memory is
56 allocated, the standard constructor does not take care of the resizing.
57 Examples for this are the STL containers like `vector<double>`.
59 The most easy way of getting your own state type to work with odeint is to use
60 a fixed size state, base calculations on the range_algebra and provide the
61 following functionality:
63 [[Name] [Expression] [Type] [Semantics]]
64 [[Construct State] [`State x()`] [`void`] [Creates an instance of `State`
65 and allocates memory.] ]
66 [[Begin of the sequence] [boost::begin(x)] [Iterator] [Returns an iterator
67 pointing to the begin of the sequence]]
68 [[End of the sequence] [boost::end(x)] [Iterator] [Returns an iterator
69 pointing to the end of the sequence]]
72 [warning If your state type does not allocate memory by default construction,
73 you [*must define it as resizeable] and provide resize functionality (see
74 below). Otherwise segmentation faults will occur.]
76 So fixed sized arrays supported by __boost_range immediately work with odeint.
77 For dynamically sized arrays one has to additionally supply the resize
79 First, the state has to be tagged as resizeable by specializing the struct
80 `is_resizeable` which consists of one typedef and one bool value:
82 [[Name] [Expression] [Type] [Semantics]]
83 [[Resizability] [`is_resizeable<State>::type`]
84 [`boost::true_type` or `boost::false_type`]
85 [Determines resizeability of the state type, returns `boost::true_type` if
86 the state is resizeable.]]
87 [[Resizability] [`is_resizeable<State>::value`]
89 [Same as above, but with `bool` value.]]
92 Defining `type` to be `true_type` and `value` as `true` tells odeint that your
94 By default, odeint now expects the support of `boost::size(x)` and a
95 `x.resize( boost::size(y) )` member function for resizing:
97 [[Name] [Expression] [Type] [Semantics]]
98 [[Get size] [`boost::size( x )`]
99 [`size_type`] [Returns the current size of x.]]
100 [[Resize] [`x.resize( boost::size( y ) )`]
101 [`void`] [Resizes x to have the same size as y.]]
104 [section Using the container interface]
105 [import ../examples/my_vector.cpp]
107 As a first example we take the most simple case and implement our own vector
108 `my_vector` which will provide a container interface.
109 This makes __boost_range working out-of-box.
110 We add a little functionality to our vector which makes it allocate some
111 default capacity by construction.
112 This is helpful when using resizing as then a resize can be assured to not
113 require a new allocation.
118 The only thing that has to be done other than defining is thus declaring
119 my_vector as resizeable:
121 [my_vector_resizeable]
123 If we wouldn't specialize the `is_resizeable` template, the code would still
124 compile but odeint would not adjust the size of temporary internal instances
125 of my_vector and hence try to fill zero-sized vectors resulting in
127 The full example can be found in [github_link examples/my_vector.cpp my_vector.cpp]
133 If your state type does work with __boost_range, but handles resizing
134 differently you are required to specialize two implementations used by odeint
135 to check a state's size and to resize:
137 [[Name] [Expression] [Type] [Semantics]]
138 [[Check size] [`same_size_impl<State,State>::same_size(x , y)`]
139 [`bool`] [Returns true if the size of x equals the size of y.]]
140 [[Resize] [`resize_impl<State,State>::resize(x , y)`]
141 [`void`] [Resizes x to have the same size as y.]]
144 As an example we will use a `std::list` as state type in odeint.
145 Because `std::list` is not supported by `boost::size` we have to replace the
146 same_size and resize implementation to get list to work with odeint.
147 The following code shows the required template specializations:
149 [import ../examples/list_lattice.cpp]
153 With these definitions odeint knows how to resize `std::list`s and so they can
154 be used as state types.
155 A complete example can be found in [github_link examples/list_lattice.cpp list_lattice.cpp].
161 [section Algebras and Operations]
163 To provide maximum flexibility odeint is implemented in a highly modularized
164 way. This means it is possible to change the underlying mathematical
165 operations without touching the integration algorithms.
166 The fundamental mathematical operations are those of a vector space, that is
167 addition of `state_types` and multiplication of `state_type`s with a scalar
168 (`time_type`). In odeint this is realized in two concepts: _Algebra_ and
170 The standard way how this works is by the range algebra which provides
171 functions that apply a specific operation to each of the individual elements
172 of a container based on the __boost_range library.
173 If your state type is not supported by __boost_range there are several
174 possibilities to tell odeint how to do algebraic operations:
176 * Implement `boost::begin` and `boost::end` for your state type so it works
178 * Implement vector-vector addition operator `+` and scalar-vector
179 multiplication operator `*` and use the non-standard `vector_space_algebra`.
180 * Implement your own algebra that implements the required functions.
184 In the following example we will try to use the `gsl_vector` type from __gsl (GNU
185 Scientific Library) as state type in odeint.
186 We will realize this by implementing a wrapper around the gsl_vector that
187 takes care of construction/destruction.
188 Also, __boost_range is extended such that it works with `gsl_vector`s as well
189 which required also the implementation of a new `gsl_iterator`.
191 [note odeint already includes all the code presented here, see [github_link
192 boost/numeric/odeint/external/gsl/gsl_wrapper.hpp gsl_wrapper.hpp], so `gsl_vector`s
193 can be used straight out-of-box.
194 The following description is just for educational purpose.]
196 The GSL is a C library, so `gsl_vector` has neither constructor, nor
197 destructor or any `begin` or `end` function, no iterators at all.
198 So to make it work with odeint plenty of things have to be implemented.
199 Note that all of the work shown here is already included in odeint, so using
200 `gsl_vector`s in odeint doesn't require any further adjustments.
201 We present it here just as an educational example.
202 We start with defining appropriate constructors and destructors.
203 This is done by specializing the `state_wrapper` for `gsl_vector`.
204 State wrappers are used by the steppers internally to create and manage
205 temporary instances of state types:
209 struct state_wrapper< gsl_vector* >
211 typedef double value_type;
212 typedef gsl_vector* state_type;
213 typedef state_wrapper< gsl_vector* > state_wrapper_type;
219 m_v = gsl_vector_alloc( 1 );
222 state_wrapper( const state_wrapper_type &x )
224 resize( m_v , x.m_v );
225 gsl_vector_memcpy( m_v , x.m_v );
231 gsl_vector_free( m_v );
237 This `state_wrapper` specialization tells odeint how gsl_vectors are created,
238 copied and destroyed.
239 Next we need resizing, this is required because gsl_vectors are dynamically
243 struct is_resizeable< gsl_vector* >
245 typedef boost::true_type type;
246 const static bool value = type::value;
250 struct same_size_impl< gsl_vector* , gsl_vector* >
252 static bool same_size( const gsl_vector* x , const gsl_vector* y )
254 return x->size == y->size;
259 struct resize_impl< gsl_vector* , gsl_vector* >
261 static void resize( gsl_vector* x , const gsl_vector* y )
263 gsl_vector_free( x );
264 x = gsl_vector_alloc( y->size );
269 Up to now, we defined creation/destruction and resizing, but gsl_vectors also
270 don't support iterators, so we first implement a gsl iterator:
274 * defines an iterator for gsl_vector
276 class gsl_vector_iterator
277 : public boost::iterator_facade< gsl_vector_iterator , double ,
278 boost::random_access_traversal_tag >
282 gsl_vector_iterator( void ): m_p(0) , m_stride( 0 ) { }
283 explicit gsl_vector_iterator( gsl_vector *p ) : m_p( p->data ) , m_stride( p->stride ) { }
284 friend gsl_vector_iterator end_iterator( gsl_vector * );
288 friend class boost::iterator_core_access;
289 friend class const_gsl_vector_iterator;
291 void increment( void ) { m_p += m_stride; }
292 void decrement( void ) { m_p -= m_stride; }
293 void advance( ptrdiff_t n ) { m_p += n*m_stride; }
294 bool equal( const gsl_vector_iterator &other ) const { return this->m_p == other.m_p; }
295 bool equal( const const_gsl_vector_iterator &other ) const;
296 double& dereference( void ) const { return *m_p; }
302 A similar class exists for the `const` version of the iterator.
303 Then we have a function returning the end iterator (similarly for `const` again):
305 gsl_vector_iterator end_iterator( gsl_vector *x )
307 gsl_vector_iterator iter( x );
308 iter.m_p += iter.m_stride * x->size;
313 Finally, the bindings for __boost_range are added:
316 inline gsl_vector_iterator range_begin( gsl_vector *x )
318 return gsl_vector_iterator( x );
322 inline gsl_vector_iterator range_end( gsl_vector *x )
324 return end_iterator( x );
327 Again with similar definitions for the `const` versions.
328 This eventually makes odeint work with gsl vectors as state types.
329 The full code for these bindings is found in [github_link
330 boost/numeric/odeint/external/gsl/gsl_wrapper.hpp gsl_wrapper.hpp].
331 It might look rather complicated but keep in mind that gsl is a pre-compiled C
336 [section Vector Space Algebra]
338 As seen above, the standard way of performing algebraic operations on
339 container-like state types in odeint is to iterate through the elements of the
340 container and perform the operations element-wise on the underlying value type.
341 This is realized by means of the `range_algebra` that uses __boost_range for
342 obtaining iterators of the state types.
343 However, there are other ways to implement the algebraic operations on
344 containers, one of which is defining the addition/multiplication operators for
345 the containers directly and then using the `vector_space_algebra`.
346 If you use this algebra, the following operators have to be defined for the
350 [[Name] [Expression] [Type] [Semantics]]
351 [[Addition] [`x + y`] [`state_type`] [Calculates the vector sum 'x+y'.]]
352 [[Assign addition] [`x += y`] [`state_type`] [Performs x+y in place.]]
353 [[Scalar multiplication] [`a * x `] [`state_type`] [Performs multiplication of vector x with scalar a.]]
354 [[Assign scalar multiplication] [`x *= a`] [`state_type`] [Performs in-place multiplication of vector x with scalar a.]]
357 Defining these operators makes your state type work with any basic Runge-Kutta
359 However, if you want to use step-size control, some more functionality is
361 Specifically, operations like
362 [' max[sub i]( |err[sub i]| / (alpha * |s[sub i]|) )]
363 have to be performed.
364 ['err] and ['s] are state_types, alpha is a scalar.
365 As you can see, we need element wise absolute value and division as well as an
366 reduce operation to get the maximum value.
367 So for controlled steppers the following things have to be implemented:
370 [[Name] [Expression] [Type] [Semantics]]
371 [[Division] [`x / y`] [`state_type`] [Calculates the element-wise division 'x/y']]
372 [[Absolute value] [`abs( x )`] [`state_type`] [Element wise absolute value]]
373 [[Reduce] [`vector_space_reduce_impl< state_type >::reduce( state , operation , init )`] [`value_type`]
374 [Performs the `operation` for subsequently each element of `state` and returns the aggregate value.
377 `init = operator( init , state[0] );`
379 `init = operator( init , state[1] )`
388 [section Boost.Ublas]
389 As an example for the employment of the `vector_space_algebra` we will adopt
390 `ublas::vector` from __ublas to work as a state type in odeint.
391 This is particularly easy because `ublas::vector` supports vector-vector
392 addition and scalar-vector multiplication described above as well as `boost::size`.
393 It also has a resize member function so all that has to be done in this case
394 is to declare resizability:
396 [import ../examples/ublas/lorenz_ublas.cpp]
400 Now ublas::vector can be used as state type for simple Runge-Kutta steppers
401 in odeint by specifying the `vector_space_algebra` as algebra in the template
402 parameter list of the stepper.
403 The following code shows the corresponding definitions:
407 Note again, that we haven't supported the requirements for controlled steppers,
408 but only for simple Runge-Kutta methods.
409 You can find the full example in [github_link
410 examples/ublas/lorenz_ublas.cpp lorenz_ublas.cpp].
417 [import ../examples/lorenz_point.cpp]
419 Here we show how to implement the required operators on a state type.
420 As example we define a new class `point3D` representing a three-dimensional
421 vector with components x,y,z and define addition and scalar multiplication
423 We use __boost_operators to reduce the amount of code to be written.
424 The class for the point type looks as follows:
428 By deriving from __boost_operators classes we don't have to define outer class
429 operators like `operator+( point3D , point3D )` because that is taken care of
430 by the operators library.
431 Note that for simple Runge-Kutta schemes (like `runge_kutta4`) only the `+`
432 and `*` operators are required.
433 If, however, a controlled stepper is used one also needs to specify the
434 division operator `/` because calculation of the error term involves an
435 element wise division of the state types.
436 Additionally, controlled steppers require an `abs` function calculating the
437 element-wise absolute value for the state type:
441 Finally, we have to provide a specialization to calculate the infintity norm of a state:
445 Again, note that the two last steps were only required if you want to use
447 For simple steppers definition of the simple `+=` and `*=` operators are
449 Having defined such a point type, we can easily perform the integration on a Lorenz
450 system by explicitely configuring the `vector_space_algebra` in the stepper's
451 template argument list:
455 The whole example can be found in [github_link
456 examples/lorenz_point.cpp lorenz_point.cpp]
458 [note For the most `state_types`, odeint is able to automatically determine
459 the correct algebra and operations. But if you want to use your own `state_type`, as in this
460 example with `point3D`, you have to manually configure the right
461 algebra/operations, unless your `state_type` works with the default choice of
462 `range_algebra` and `default_operations`.]
468 gsl_vector, gsl_matrix, ublas::matrix, blitz::matrix, thrust
470 [section Adapt your own operations]