ceph/src/boost/libs/numeric/odeint/doc/tutorial_special_topics.qbk

   1 [/============================================================================
   2   Boost.odeint
   3
   4   Copyright 2011-2012 Karsten Ahnert
   5   Copyright 2011-2013 Mario Mulansky
   6   Copyright 2012 Sylwester Arabas
   7
   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 =============================================================================/]
  12
  13
  14 [/ [section Special topics] /]
  15
  16 [section Complex state types]
  17
  18 [import ../examples/stuart_landau.cpp]
  19
  20 Thus far we have seen several examples defined for real values.
  21 odeint can handle complex state types, hence ODEs which are defined on complex
  22 vector spaces, as well. An example is the Stuart-Landau oscillator
  23
  24 ['d __Psi / dt = ( 1 + i __eta ) __Psi + ( 1 + i __alpha ) | __Psi |[super 2] __Psi ]
  25
  26 where ['__Psi] and ['i] is a complex variable. The definition of this ODE in C++
  27 using complex< double > as a state type may look as follows
  28
  29 [stuart_landau_system_function]
  30
  31 One can also use a function instead of a functor to implement it
  32
  33 [stuart_landau_system_function_alternative]
  34
  35 We strongly recommend to use the first ansatz. In this case you have explicit control over the parameters of the system and are not restricted to use global variables to parametrize the oscillator.
  36
  37 When choosing the stepper type one has to account for the "unusual" state type:
  38 it is a single `complex<double>` opposed to the vector types used in the
  39 previous examples. This means that no iterations over vector elements have to
  40 be performed inside the stepper algorithm. Odeint already detects that and
  41 automatically uses the `vector_space_algebra` for computation.
  42 You can enforce this by supplying additional template arguments to the stepper
  43 including the `vector_space_algebra`. Details on the usage of algebras can be
  44 found in the section __adapt_state_types.
  45
  46 [stuart_landau_integration]
  47
  48 The full cpp file for the Stuart-Landau example can be found here [github_link
  49 examples/stuart_landau.cpp stuart_landau.cpp]
  50
  51 [endsect]
  52
  53 [section Lattice systems]
  54
  55 [import ../examples/fpu.cpp]
  56
  57
  58 odeint can also be used to solve ordinary differential equations defined on lattices. A prominent example is the Fermi-Pasta-Ulam system __fpu_scholarpedia_ref. It is a Hamiltonian system of nonlinear coupled harmonic oscillators. The Hamiltonian is
  59
  60 [' H = __Sigma[subl i] p[subl i][super 2]/2 + 1/2 ( q[subl i+1] - q[subl i] )^2 + __beta / 4 ( q[subl i+1] - q[subl i] )^4 ]
  61
  62 Remarkably, the Fermi-Pasta-Ulam system was the first numerical experiment to be implemented on a computer. It was studied at Los Alamos in 1953 on one of the first computers (a MANIAC I) and it triggered a whole new tree of mathematical and physical science.
  63
  64 Like the __tut_solar_system, the FPU is solved again by a symplectic solver, but in this case we can speed up the computation because the ['q] components trivially reduce to ['dq[subl i] / dt = p[subl i]]. odeint is capable of doing this performance improvement. All you have to do is to call the symplectic solver with an state function for the ['p] components. Here is how this function looks like
  65
  66 [fpu_system_function]
  67
  68 You can also use `boost::array< double , N >` for the state type.
  69
  70 Now, you have to define your initial values and perform the integration:
  71
  72 [fpu_integration]
  73
  74 The observer uses a reference to the system object to calculate the local energies:
  75
  76 [fpu_observer]
  77
  78 The full cpp file for this FPU example can be found here [github_link examples/fpu.cpp fpu.cpp]
  79
  80 [endsect]
  81
  82 [section Ensembles of oscillators]
  83
  84 [import ../examples/phase_oscillator_ensemble.cpp]
  85
  86 Another important high dimensional system of coupled ordinary differential equations is an ensemble of ['N] all-to-all coupled phase oscillators __synchronization_pikovsky_ref. It is defined as
  87
  88 [' d__phi[subl k] / dt = __omega[subl k] + __epsilon / N __Sigma[subl j] sin( __phi[subl j] - __phi[subl k] )]
  89
  90 The natural frequencies ['__omega[subl i]] of each oscillator follow some distribution and ['__epsilon] is the coupling strength. We choose here a Lorentzian distribution for ['__omega[subl i]]. Interestingly a phase transition can be observed if the coupling strength exceeds a critical value. Above this value synchronization sets in and some of the oscillators oscillate with the same frequency despite their different natural frequencies. The transition is also called Kuramoto transition. Its behavior can be analyzed by employing the mean field of the phase
  91
  92 ['Z = K e[super i __Theta] = 1 / N __Sigma[subl k]e[super i __phi[subl k]]]
  93
  94 The definition of the system function is now a bit more complex since we also need to store the individual frequencies of each oscillator.
  95
  96 [phase_oscillator_ensemble_system_function]
  97
  98 Note, that we have used ['Z] to simplify the equations of motion. Next, we create an observer which computes the value of ['Z] and we record ['Z] for different values of ['__epsilon].
  99
 100 [phase_oscillator_ensemble_observer]
 101
 102 Now, we do several integrations for different values of ['__epsilon] and record ['Z]. The result nicely confirms the analytical result of the phase transition, i.e. in our example the standard deviation of the Lorentzian is 1 such that the transition will be observed at ['__epsilon = 2].
 103
 104 [phase_oscillator_ensemble_integration]
 105
 106 The full cpp file for this example can be found here [github_link examples/phase_oscillator_ensemble.cpp phase_oscillator_ensemble.cpp]
 107
 108 [endsect]
 109
 110 [section Using boost::units]
 111
 112 [import ../examples/harmonic_oscillator_units.cpp]
 113
 114 odeint also works well with __boost_units - a library for compile time unit
 115 and dimension analysis. It works by decoding unit information into the types
 116 of values. For a one-dimensional unit you can just use the Boost.Unit types as
 117 state type, deriv type and time type and hand the `vector_space_algebra` to
 118 the stepper definition and everything works just fine:
 119
 120 ```
 121 typedef units::quantity< si::time , double > time_type;
 122 typedef units::quantity< si::length , double > length_type;
 123 typedef units::quantity< si::velocity , double > velocity_type;
 124
 125 typedef runge_kutta4< length_type , double , velocity_type , time_type ,
 126                       vector_space_algebra > stepper_type;
 127 ```
 128
 129 If you want to solve more-dimensional problems the individual entries
 130 typically have different units. That means that the `state_type` is now
 131 possibly heterogeneous, meaning that every entry might have a different type.
 132 To solve this problem, compile-time sequences from __boost_fusion can be used.
 133
 134 To illustrate how odeint works with __boost_units we use the harmonic oscillator as primary example. We start with defining all quantities
 135
 136 [units_define_basic_quantities]
 137
 138 Note, that the `state_type` and the `deriv_type` are now a compile-time fusion
 139 sequences. `deriv_type` represents ['x'] and is now different from the state
 140 type as it has different unit definitions. Next, we define the ordinary
 141 differential equation which is completely equivalent to the example in __tut_harmonic_oscillator:
 142
 143 [units_define_ode]
 144
 145 Next, we instantiate an appropriate stepper. We must explicitly parametrize
 146 the stepper with the `state_type`, `deriv_type`, `time_type`.
 147
 148 [units_define_stepper]
 149
 150 [note When using compile-time sequences, the iteration over vector elements is
 151 done by the `fusion_algebra`, which is automatically chosen by odeint. For
 152 more on the state types / algebras see chapter __adapt_state_types.]
 153
 154 It is quite easy but the compilation time might take very long. Furthermore, the observer is defined a bit different
 155
 156 [units_observer]
 157
 158 [caution Using __boost_units works nicely but compilation can be very time and
 159 memory consuming. For example the unit test for the usage of __boost_units in odeint take up to 4 GB
 160 of memory at compilation.]
 161
 162 The full cpp file for this example can be found here [github_link examples/harmonic_oscillator_units.cpp harmonic_oscillator_units.cpp].
 163
 164 [endsect]
 165
 166 [section Using matrices as state types]
 167
 168 [import ../examples/two_dimensional_phase_lattice.cpp]
 169
 170 odeint works well with a variety of different state types. It is not restricted to pure vector-wise types, like `vector< double >`, `array< double , N >`, `fusion::vector< double , double >`, etc. but also works with types having a different topology then simple vectors. Here, we show how odeint can be used with matrices as states type, in the next section we will show how can be used to solve ODEs defined on complex networks.
 171
 172 By default, odeint can be used with `ublas::matrix< T >` as state type for matrices. A simple example is a two-dimensional lattice of coupled phase oscillators. Other matrix types like `mtl::dense_matrix` or blitz arrays and matrices can used as well but need some kind of activation in order to work with odeint. This activation is described in following sections,
 173
 174 The definition of the system is
 175
 176 [two_dimensional_phase_lattice_definition]
 177
 178 In principle this is all. Please note, that the above code is far from being optimal. Better performance can be achieved if every interaction is only calculated once and iterators for columns and rows are used. Below are some visualizations of the evolution of this lattice equation.
 179
 180 [$phase_lattice_2d_0000.jpg] [$phase_lattice_2d_0100.jpg] [$phase_lattice_2d_1000.jpg]
 181
 182 The full cpp for this example can be found here [github_link examples/two_dimensional_phase_lattice.cpp two_dimensional_phase_lattice.cpp].
 183
 184 [endsect]
 185
 186 [/
 187 [section Partial differential equations]
 188 To be continued:
 189 *Wave equation
 190 *KdV
 191 *Ginzburg-Landau
 192 [endsect]
 193 [section Ordinary differential equations on networks]
 194 to be continued
 195 [endsect]
 196 ]
 197
 198 [section Using arbitrary precision floating point types]
 199
 200 [import ../examples/multiprecision/lorenz_mp.cpp]
 201
 202 Sometimes one needs results with higher precision than provided by the
 203 standard floating point types.
 204 As odeint allows to configure the fundamental numerical type, it is well
 205 suited to be run with arbitrary precision types.
 206 Therefore, one only needs a library that provides a type representing values
 207 with arbitrary precision and the fundamental operations for those values.
 208 __boost_multiprecision is a boost library that does exactly this.
 209 Making use of __boost_multiprecision to solve odes with odeint is very simple,
 210 as the following example shows.
 211
 212 Here we use `cpp_dec_float_50` as the fundamental value type, which ensures
 213 exact computations up to 50 decimal digits.
 214
 215
 216 [mp_lorenz_defs]
 217
 218 As exemplary ODE again the lorenz system is chosen, but here we have to make
 219 sure all constants are initialized as high precision values.
 220
 221 [mp_lorenz_rhs]
 222
 223 The actual integration then is straight forward:
 224
 225 [mp_lorenz_int]
 226
 227 The full example can be found at [github_link examples/multiprecision/lorenz_mp.cpp lorenz_mp.cpp].
 228 Another example that compares the accuracy of the high precision type with
 229 standard double can be found at [github_link examples/multiprecision/cmp_precision.cpp cmp_precision.cpp].
 230
 231 Furthermore, odeint can also be run with other multiprecision libraries,
 232 e.g. [@http://gmplib.org/ gmp].
 233 An example for this is given in [github_link examples/gmpxx/lorenz_gmpxx.cpp lorenz_gmpxx.cpp].
 234
 235 [endsect]
 236
 237 [section Self expanding lattices]
 238
 239 [import ../examples/resizing_lattice.cpp]
 240
 241 odeint supports changes of the state size during integration if a state_type
 242 is used which can be resized, like `std::vector`.
 243 The adjustment of the state's size has to be done from outside and the stepper
 244 has to be instantiated with `always_resizer` as the template argument for the
 245 `resizer_type`.
 246 In this configuration, the stepper checks for changes in the state size and
 247 adjust it's internal storage accordingly.
 248
 249 We show this for a Hamiltonian system of nonlinear, disordered oscillators with nonlinear nearest neighbor coupling.
 250
 251 The system function is implemented in terms of a class that also provides functions for calculating the energy.
 252 Note, that this class stores the random potential internally which is not resized, but rather a start index is kept which should be changed whenever the states' size change.
 253
 254 [resizing_lattice_system_class]
 255
 256 The total size we allow is 1024 and we start with an initial state size of 60.
 257
 258 [resizing_lattice_initialize]
 259
 260 The lattice gets resized whenever the energy distribution comes close to the borders `distr[10] > 1E-150`, `distr[distr.size()-10] > 1E-150`.
 261 If we increase to the left, `q` and `p` have to be rotated because their resize function always appends at the end.
 262 Additionally, the start index of the potential changes in this case.
 263
 264 [resizing_lattice_steps_loop]
 265
 266 The `do_resize` function simply calls `vector.resize` of `q` , `p` and `distr`.
 267
 268 [resizing_lattice_resize_function]
 269
 270 The full example can be found in [github_link examples/resizing_lattice.cpp resizing_lattice.cpp]
 271
 272 [endsect]
 273
 274 [/ [endsect] /]