1[/============================================================================ 2 Boost.odeint 3 4 Copyright 2011-2013 Karsten Ahnert 5 Copyright 2011-2012 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 15[section Getting started] 16 17[section Overview] 18 19odeint is a library for solving initial value problems (IVP) of ordinary 20differential equations. Mathematically, these problems are formulated as 21follows: 22 23['x'(t) = f(x,t)], ['x(0) = x0]. 24 25['x] and ['f] can be vectors and the solution is some function ['x(t)] fulfilling both equations above. In the following we will refer to ['x'(t)] also `dxdt` which is also our notation for the derivative in the source code. 26 27Ordinary differential equations occur nearly everywhere in natural sciences. For example, the whole Newtonian mechanics are described by second order differential equations. Be sure, you will find them in every discipline. They also occur if partial differential equations (PDEs) are discretized. Then, a system of coupled ordinary differential occurs, sometimes also referred as lattices ODEs. 28 29Numerical approximations for the solution ['x(t)] are calculated iteratively. The easiest algorithm is the Euler scheme, where starting at ['x(0)] one finds ['x(dt) = x(0) + dt f(x(0),0)]. Now one can use ['x(dt)] and obtain ['x(2dt)] in a similar way and so on. The Euler method is of order 1, that means the error at each step is ['~ dt[super 2]]. This is, of course, not very satisfying, which is why the Euler method is rarely used for real life problems and serves just as illustrative example. 30 31The main focus of odeint is to provide numerical methods implemented in a way where the algorithm is completely independent on the data structure used to represent the state /x/. 32In doing so, odeint is applicable for a broad variety of situations and it can be used with many other libraries. Besides the usual case where the state is defined as a `std::vector` or a `boost::array`, we provide native support for the following libraries: 33 34* __ublas 35* __thrust, making odeint naturally running on CUDA devices 36* gsl_vector for compatibility with the many numerical function in the GSL 37* __boost_range 38* __boost_fusion (the state type can be a fusion vector) 39* __boost_units 40* __intel_mkl for maximum performance 41* __vexcl for OpenCL 42* __boost_graph (still experimentally) 43 44In odeint, the following algorithms are implemented: 45 46[include stepper_table.qbk] 47 48[endsect] 49 50 51 52[section Usage, Compilation, Headers] 53 54odeint is a header-only library, no linking against pre-compiled code is required. It can be included by 55 56`` 57#include <boost/numeric/odeint.hpp> 58`` 59 60which includes all headers of the library. All functions and classes from odeint live in the namespace 61 62`` 63 using namespace boost::numeric::odeint; 64`` 65 66It is also possible to include only parts of the library. This is the recommended way since it saves a lot of compilation time. 67 68* `#include <boost/numeric/odeint/stepper/XYZ.hpp>` - the include path for all steppers, XYZ is a placeholder for a stepper. 69* `#include <boost/numeric/odeint/algebra/XYZ.hpp>` - all algebras. 70* `#include <boost/numeric/odeint/util/XYZ.hpp>` - the utility functions like `is_resizeable`, `same_size`, or `resize`. 71* `#include <boost/numeric/odeint/integrate/XYZ.hpp>` - the integrate routines. 72* `#include <boost/numeric/odeint/iterator/XYZ.hpp>` - the range and iterator functions. 73* `#include <boost/numeric/odeint/external/XYZ.hpp>` - any binders to external libraries. 74 75 76 77[endsect] 78 79[section Short Example] 80 81Imaging, you want to numerically integrate a harmonic oscillator with 82friction. The equations of motion are given by ['x'' = -x + __gamma x']. 83Odeint only deals with first order ODEs that have no higher derivatives than 84x' involved. However, any higher order ODE can be transformed to a system of 85first order ODEs by introducing the new variables ['q=x] and ['p=x'] such that ['w=(q,p)]. To apply numerical integration one first has to design the right hand side of the equation ['w' = f(w) = (p,-q+__gamma p)]: 86 87[import ../examples/harmonic_oscillator.cpp] 88[rhs_function] 89 90Here we chose `vector<double>` as the state type, but others are also 91possible, for example `boost::array<double,2>`. odeint is designed in such a 92way that you can easily use your own state types. Next, the ODE is defined 93which is in this case a simple function calculating ['f(x)]. The parameter 94signature of this function is crucial: the integration methods will always 95call them in the form `f(x, dxdt, t)` (there are exceptions for some special routines). So, even if there is no explicit time dependence, one has to define `t` as a function parameter. 96 97Now, we have to define the initial state from which the integration should start: 98 99[state_initialization] 100 101For the integration itself we'll use the `integrate` function, which is a convenient way to get quick results. It is based on the error-controlled `runge_kutta54_cash_karp` stepper (5th order) and uses adaptive step-size. 102 103[integration] 104 105The integrate function expects as parameters the rhs of the ode as defined 106above, the initial state `x`, the start-and end-time of the integration as 107well as the initial time step=size. Note, that `integrate` uses an adaptive step-size during 108the integration steps so the time points will not be equally spaced. The 109integration returns the number of steps that were applied and updates x which 110is set to the approximate solution of the ODE at the end of integration. 111 112It is also possible to represent the ode system as a class. The 113rhs must then be implemented as a functor - a class with an overloaded function call operator: 114 115[rhs_class] 116 117which can be used via 118 119[integration_class] 120 121In order to observe the solution 122during the integration steps all you have to do is 123to provide a reasonable observer. An example is 124 125[integrate_observer] 126 127which stores the intermediate steps in a container. 128Note, the argument structure of the ()-operator: odeint calls the observer 129exactly in this way, providing the current state and time. Now, you only have 130to pass this container to the integration function: 131 132[integrate_observ] 133 134That is all. You can use functional libraries like __boost_lambda or __boost_phoenix to ease the creation of observer functions. 135 136The full cpp file for this example can be found here: [github_link examples/harmonic_oscillator.cpp harmonic_oscillator.cpp] 137 138 139 140[endsect] 141 142[endsect] 143