diff --git a/doc/boost_numeric_odeint/concepts.html b/doc/boost_numeric_odeint/concepts.html deleted file mode 100644 index 6c6e772c..00000000 --- a/doc/boost_numeric_odeint/concepts.html +++ /dev/null @@ -1,54 +0,0 @@ - - - -Concepts - - - - - - - - -
-PrevUpHomeNext -
-
-

-Concepts -

-
-
System
-
Symplectic - System
-
Simple - Symplectic System
-
Implicit - System
-
Stepper
-
Error Stepper
-
Controlled - Stepper
-
Dense - Output Stepper
-
State - Algebra Operations
-
State Wrapper
-
- -
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/controlled_stepper.html b/doc/boost_numeric_odeint/concepts/controlled_stepper.html deleted file mode 100644 index 7e317d74..00000000 --- a/doc/boost_numeric_odeint/concepts/controlled_stepper.html +++ /dev/null @@ -1,227 +0,0 @@ - - - -Controlled Stepper - - - - - - - - -
-PrevUpHomeNext -
-
-

-Controlled - Stepper -

-

- This concept specifies the interface a controlled stepper has to fulfill - to be used within integrate - functions. -

-
- - Description -
-

- A controlled stepper following this Controlled Stepper concept provides the - possibility to perform one step of the solution x(t) - of an ODE with step-size dt to obtain x(t+dt) - with a given step-size dt. Depending on an error estimate - of the solution the step might be rejected and a smaller step-size is suggested. -

-
- - Associated - types -
-
-
- - Notation -
-
-

-
-
ControlledStepper
-

- A type that is a model of Controlled Stepper -

-
State
-

- A type representing the state x of the ODE -

-
Time
-

- A type representing the time t of the ODE -

-
stepper
-

- An object of type ControlledStepper -

-
x
-

- Object of type State -

-
t, dt
-

- Objects of type Time -

-
sys
-

- An object defining the ODE, should be a model of System, - Symplectic - System, Simple - Symplectic System or Implicit - System. -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Do step -

- 		-

-

-
stepper.try_step( sys , x , t , dt )
-

-

- 		-

- controlled_step_result -

- 		-

- Tries one step of step size dt. - If the step was successful, success - is returned, the resulting state is written to x, - the new time is stored in t - and dt now contains - a new (possibly larger) step-size for the next step. If the error - was too big, rejected - is returned and the results are neglected - x - and t are unchanged - and dt now contains - a reduced step-size to be used for the next try. -

-
-
- - Models -
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/dense_output_stepper.html b/doc/boost_numeric_odeint/concepts/dense_output_stepper.html deleted file mode 100644 index b7039d58..00000000 --- a/doc/boost_numeric_odeint/concepts/dense_output_stepper.html +++ /dev/null @@ -1,311 +0,0 @@ - - - -Dense Output Stepper - - - - - - - - -
-PrevUpHomeNext -
-
-

-Dense - Output Stepper -

-

- This concept specifies the interface a dense output stepper has to fulfill - to be used within integrate - functions. -

-
- - Description -
-

- A dense output stepper following this Dense Output Stepper concept provides - the possibility to perform a single step of the solution x(t) - of an ODE to obtain x(t+dt). The step-size dt might be adjusted automatically due - to error control. Dense output steppers also can interpolate the solution - to calculate the state x(t') at any point t - <= t' <= t+dt. -

-
- - Associated - types -
-
-
- - Notation -
-
-

-
-
Stepper
-

- A type that is a model of Dense Output Stepper -

-
State
-

- A type representing the state x of the ODE -

-
stepper
-

- An object of type Stepper -

-
x0, x
-

- Object of type State -

-
t0, dt0, t
-

- Objects of type Stepper::time_type -

-
sys
-

- An object defining the ODE, should be a model of System, - Symplectic - System, Simple - Symplectic System or Implicit - System. -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Initialize integration -

- 		-

- stepper.initialize( - x0 , - t0 , - dt0 ) -

- 		-

- void -

- 		-

- Initializes the stepper with initial values x0, - t0 and dt0. -

-
-

- Do step -

- 		-

- stepper.do_step( - sys ) -

- 		-

- std::pair< - Stepper::time_type , - Stepper::time_type > -

- 		-

- Performs one step using the ODE defined by sys. - The step-size might be changed internally due to error control. - This function returns a pair containing t - and t+dt representing the interval - for which interpolation can be performed. -

-
-

- Do interpolation -

- 		-

- stepper.calc_state( - t_inter , - x ) -

- 		-

- void -

- 		-

- Performs the interpolation to calculate /x(tinter/) where /t <= - tinter <= t+dt/. -

-
-

- Get current time -

- 		-

- stepper.current_time() -

- 		-

- const Stepper::time_type& -

- 		-

- Returns the current time t+dt of the stepper, - that is the end time of the last step and the starting time for - the next call of do_step -

-
-

- Get current state -

- 		-

- stepper.current_state() -

- 		-

- const Stepper::state_type& -

- 		-

- Returns the current state of the stepper, that is x(t+dt), - the state at the time returned by stepper.current_time() -

-
-
- - Models -
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/error_stepper.html b/doc/boost_numeric_odeint/concepts/error_stepper.html deleted file mode 100644 index bf866148..00000000 --- a/doc/boost_numeric_odeint/concepts/error_stepper.html +++ /dev/null @@ -1,361 +0,0 @@ - - - -Error Stepper - - - - - - - - -
-PrevUpHomeNext -
-
-

-Error Stepper -

-

- This concepts specifies the interface an error stepper has to fulfill to - be used within a ControlledErrorStepper. An error stepper must always fullfil - the stepper concept. This can trivially implemented by -

-

-

-
template< class System >
-error_stepper::do_step( System sys , state_type &x , time_type t , time_type dt )
-{
-    state_type xerr;
-    // allocate xerr
-    do_step( sys , x , t , dt , xerr );
-}
-
-

-

-
- - Description -
-

- An error stepper following this Error Stepper concept is capable of doing - one step of the solution x(t) of an ODE with step-size - dt to obtain x(t+dt) and also computing - an error estimate xerr of the result. Error Steppers - can be Runge Kutta steppers, symplectic steppers as well as implicit steppers. - Based on the stepper type, the ODE is defined as System, - Symplectic - System, Simple - Symplectic System or Implicit - System. -

-
- - Refinement - of -
-
-
- - Associated - types -
-
-
- - Notation -
-
-

-
-
ErrorStepper
-

- A type that is a model of Error Stepper -

-
State
-

- A type representing the state x of the ODE -

-
Error
-

- A type representing the error calculated by the stepper, usually same - as State -

-
Time
-

- A type representing the time t of the ODE -

-
stepper
-

- An object of type ErrorStepper -

-
x
-

- Object of type State -

-
xerr
-

- Object of type Error -

-
t, dt
-

- Objects of type Time -

-
sys
-

- An object defining the ODE, should be a model of either System, - Symplectic - System, Simple - Symplectic System or Implicit - System. -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Get the stepper order -

- 		-

- stepper.order() -

- 		-

- order_type -

- 		-

- Returns the order of the stepper for one step without error estimation. -

-
-

- Get the stepper order -

- 		-

- stepper.stepper_order() -

- 		-

- order_type -

- 		-

- Returns the order of the stepper for one error estimation step - which is used for error calculation. -

-
-

- Get the error order -

- 		-

- stepper.errorr_order() -

- 		-

- order_type -

- 		-

- Returns the order of the error step which is used for error calculation. -

-
-

- Do step -

- 		-

- stepper.do_step( - sys , - x , - t , - dt ) -

- 		-

- void -

- 		-

- Performs one step of step size dt. - The newly obtained state is written in-place to x. -

-
-

- Do step with error estimation -

- 		-

- stepper.do_step( - sys , - x , - t , - dt , - xerr ) -

- 		-

- void -

- 		-

- Performs one step of step size dt - with error estimation. The newly obtained state is written in-place - to x and the estimated - error to xerr. -

-
-
- - Models -
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/implicit_system.html b/doc/boost_numeric_odeint/concepts/implicit_system.html deleted file mode 100644 index 51be7ef5..00000000 --- a/doc/boost_numeric_odeint/concepts/implicit_system.html +++ /dev/null @@ -1,175 +0,0 @@ - - - -Implicit System - - - - - - - - -
-PrevUpHomeNext -
-
-

-Implicit - System -

-
- - Description -
-

- This concept describes how to define a ODE that can be solved by an implicit - routine. Implicit routines need not only the function f(x,t) - but also the Jacobian df/dx = A(x,t). A - is a matrix and implicit routines need to solve the linear problem Ax - = b. In odeint this is implemented with use of Boost.UBlas, - therefore, the state_type implicit routines is ublas::vector - and the matrix is defined as ublas::matrix. -

-
- - Notation -
-
-

-
-
System
-

- A type that is a model of Implicit_System -

-
Time
-

- A type representing the time of the ODE -

-
sys
-

- An object of type System -

-
x
-

- Object of type ublas::vector -

-
dxdt
-

- Object of type ublas::vector -

-
jacobi
-

- Object of type ublas::matrix -

-
t
-

- Object of type Time -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Calculate dx/dt := f(x,t) -

- 		-

- sys.first( - x , - dxdt , - t ) -

- 		-

- void -

- 		-

- Calculates f(x,t), - the result is stored into dxdt -

-
-

- Calculate A := df/dx (x,t) -

- 		-

- sys.second( - x , - jacobi , - t ) -

- 		-

- void -

- 		-

- Calculates the Jacobian of f at x,t, - the result is stored into jacobi -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/simple_symplectic_system.html b/doc/boost_numeric_odeint/concepts/simple_symplectic_system.html deleted file mode 100644 index eab3bae7..00000000 --- a/doc/boost_numeric_odeint/concepts/simple_symplectic_system.html +++ /dev/null @@ -1,179 +0,0 @@ - - - -Simple Symplectic System - - - - - - - - -
-PrevUpHomeNext -
-
-

-Simple - Symplectic System -

-
- - Description -
-

- In most Hamiltonian systems the kinetic term is a quadratic term in the momentum - Hkin = p^2 / 2m and in many cases it is possible to rescale - coordinates and set m=1 which leads to a trivial equation - of motion: -

-

- q'(t) = f(p) = p. -

-

- while for p' we still have the general form -

-

- p'(t) = g(q) -

-

- As this case is very frequent we introduced a concept where only the nontrivial - equation for p' has to be provided to the symplectic - stepper. We call this concept Simple_Symplectic_System -

-
- - Notation -
-
-

-
-
System
-

- A type that is a model of Simple_Symplectic_System -

-
Coor
-

- The type of the coordinate q -

-
MomentumDeriv
-

- The type of the derivative of momentum p' -

-
sys
-

- An object that models System -

-
q
-

- Object of type Coor -

-
dpdt
-

- Object of type MomentumDeriv -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Check for pair -

- 		-

- boost::is_pair< - System >::type -

- 		-

- boost::mpl::false_ -

- 		-

- Check if System is a pair, should be evaluated to false in this - case. -

-
-

- Calculate dp/dt = g(q) -

- 		-

- sys( - q , - dpdt ) -

- 		-

- void -

- 		-

- Calculates g(q), the result is stored into - dpdt -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/state_algebra_operations.html b/doc/boost_numeric_odeint/concepts/state_algebra_operations.html deleted file mode 100644 index 6feedb65..00000000 --- a/doc/boost_numeric_odeint/concepts/state_algebra_operations.html +++ /dev/null @@ -1,827 +0,0 @@ - - - -State Algebra Operations - - - - - - - - -
-PrevUpHomeNext -
-
-

-State - Algebra Operations -

-
-
Operations
-
Algebra
-
Pre-Defined - implementations
-
Example - expressions
-
-
- - - - - -
[Note]Note

- The following does not apply to implicit steppers like implicit_euler or - Rosenbrock 4 as there the state_type - can not be changed from ublas::vector - and no algebra/operations are used. -

-
- - Description -
-

- The State, Algebra and Operations - together define a concept describing how the mathematical vector operations - required for the stepper algorithms are performed. The typical vector operation - done within steppers is -

-

- y = Σ αi xi. -

-

- The State represents the - state variable of an ODE, usually denoted with x. Algorithmically, - the state is often realized as a vector< double > or array< double , N >, - however, the genericity of odeint enables you to basically use anything as - a state type. The algorithmic counterpart of such mathematical expressions - is divided into two parts. First, the Algebra - is used to account for the vector character of the equation. In the case - of a vector as state type - this means the Algebra is - responsible for iteration over all vector elements. Second, the Operations are used to represent the actual - operation applied to each of the vector elements. So the Algebra - iterates over all elements of the States - and calls an operation taken from the Operations - for each element. This is where State, - Algebra and Operations have to work together to make - odeint running. Please have a look at the range_algebra - and default_operations to - see an example how this is implemented. -

-

- In the following we describe how State, - Algebra and Operations are used together within the - stepper implementations. -

-
-

-Operations -

-
- - Notation -
-
-

-
-
Operations
-

- The operations type -

-
Value1, ... , - ValueN
-

- Types representing the value or time type of stepper -

-
Scale
-

- Type of the scale operation -

-
scale
-

- Object of type Scale -

-
ScaleSumN
-

- Type that represents a general scale_sum operation, N - should be replaced by a number from 1 to 14. -

-
scale_sum/N/
-

- Object of type ScaleSumN, - N should be replaced by a - number from 1 to 14. -

-
ScaleSumSwap2
-

- Type of the scale sum swap operation -

-
scale_sum_swap2
-

- Object of type ScaleSumSwap2 -

-
a1, - a2, - ...
-

- Objects of type Value1, - Value2, ... -

-
y, - x1, - x2, - ...
-

- Objects of State's - value type -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Get scale operation -

- 		-

- Operations::scale< - Value > -

- 		-

- Scale -

- 		-

- Get Scale from - Operations -

-
-

- Scale constructor -

- 		-

- Scale< - Value >( - a ) -

- 		-

- Scale -

- 		-

- Constructs a Scale - object -

-
-

- Scale operation -

- 		-

- scale( - x ) -

- 		-

- void -

- 		-

- Calculates x *= a -

-
-

- Get general scale_sum - operation -

- 		-

- Operations::scale_sum/N/< Value1 , ... , ValueN > -

- 		-

- ScaleSum/N/ -

- 		-

- Get the ScaleSum/N/ type from Operations, N - should be replaced by a number from 1 to 14. -

-
-

- scale_sum constructor -

- 		-

- ScaleSum/N/< Value1 , ... , ValueN >( a1 , ... - , aN ) -

- 		-

- ScaleSum/N/ -

- 		-

- Constructs a scale_sum - object given N parameter - values with N between - 1 and 14. -

-
-

- scale_sum operation -

- 		-

- scale_sum/N/( y , x1 , ... , xN ) -

- 		-

- void -

- 		-

- Calculates y = a1*x1 - + a2*x2 - + ... - + aN*xN. - Note that this is an N+1-ary - function call. -

-
-

- Get scale sum swap operation -

- 		-

- Operations::scale_sum_swap2< - Value1 , - Value2 > -

- 		-

- ScaleSumSwap2 -

- 		-

- Get scale sum swap from operations -

-
-

- ScaleSumSwap2 - constructor -

- 		-

- ScaleSumSwap2< Value1 - , Value2 - >( a1 - , a2 - ) -

- 		-

- ScaleSumSwap2 -

- 		-

- Constructor -

-
-

- ScaleSumSwap2 - operation -

- 		-

- scale_sum_swap2( x1 - , x2 - , x3 - ) -

- 		-

- void -

- 		-

- Calculates tmp = x1, - x1 = - a1*x2 + - a2*x3 and x2 - = tmp. -

-
-
-
-

-Algebra -

-
- - Notation -
-
-

-
-
State
-

- The state type -

-
Algebra
-

- The algebra type -

-
Operation/N/
-

- An N-ary operation type, - N should be a number from - 1 to 14. -

-
algebra
-

- Object of type Algebra -

-
operation/N/
-

- Object of type Operation/N -

-
y, - x1, - x2, - ...
-

- Objects of type State -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Vector Operation with arity 2 -

- 		-

- algebra.for_each2( - y , - x , - operation2 ) -

- 		-

- void -

- 		-

- Calls operation2( y_i - , x_i - ) for each element y_i of y - and x_i of x. -

-
-

- Vector Operation with arity 3 -

- 		-

- algebra.for_each3( - y , - x1 , - x2 , - operation3 ) -

- 		-

- void -

- 		-

- Calls operation3( y_i - , x1_i - , x2_i - ) for each element y_i of y - and x1_i of - x1 and x2_i of x2. -

-
-

- Vector Operation with arity N -

- 		-

- algebra.for_each/N/( y , x1 , ... , xN , operation/N/ - ) -

- 		-

- void -

- 		-

- Calls operation/N/( y_i , x1_i , ... , xN_i ) - for each element y_i - of y and x1_i of x1 - and so on. N should be - replaced by a number between 1 and 14. -

-
-
-
-

-Pre-Defined - implementations -

-

- As standard configuration odeint uses the range_algebra - and default_operations - which suffices most situations. However, a few more possibilities exist - either to gain better performance or to ensure interoperability with other - libraries. In the following we list the existing Algebra/Operations configurations that can be - used in the steppers. -

-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- State -

- 		-

- Algebra -

- 		-

- Operations -

- 		-

- Remarks -

-
-

- Anything supporting Boost.Range, - like std::vector, std::list, - boost::array,... based on a value_type that supports operators - +,* (typically double) -

- 		-

- range_algebra -

- 		-

- default_operations -

- 		-

- Standard implementation, applicable for most typical situations. -

-
-

- boost::array based on a value_type that supports operators - +,* -

- 		-

- array_algebra -

- 		-

- default_operations -

- 		-

- Special implementation for boost::array with better performance - than range_algebra -

-
-

- Anything that defines operators + within itself and * with scalar - (Mathematically spoken, anything that is a vector space). -

- 		-

- vector_space_algebra -

- 		-

- default_operations -

- 		-

- For the use of Controlled - Stepper, the template vector_space_reduce - has to be instantiated. -

-
-

- thrust::device_vector, thrust::host_vector -

- 		-

- thrust_algebra -

- 		-

- thrust_operations -

- 		-

- For running odeint on CUDA devices by using Thrust -

-
-

- boost::array or anything which allocates - the elements in a C-like manner -

- 		-

- vector_space_algebra -

- 		-

- mkl_operations -

- 		-

- Using the Intel - Math Kernel Library in odeint for maximum performance. - Currently, only the RK4 stepper is supported. -

-
-
-
-

-Example - expressions -

-
------ - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Vector operation -

- 		-

- algebra.for_each3( - y , - x1 , - x2 , - Operations::scale_sum2< - Value1 , - Value2 >( - a1 , - a2 ) - ) -

- 		-

- void -

- 		-

- Calculates y = a1 - x1 + a2 x2 -

-
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/state_wrapper.html b/doc/boost_numeric_odeint/concepts/state_wrapper.html deleted file mode 100644 index 6669af39..00000000 --- a/doc/boost_numeric_odeint/concepts/state_wrapper.html +++ /dev/null @@ -1,278 +0,0 @@ - - - -State Wrapper - - - - - - - -
-PrevUpHome -
-
-

-State Wrapper -

-
- - Description -
-

- The State Wrapper - concept describes the way odeint creates temporary state objects to store - intermediate results within the stepper's do_step - methods. -

-
- - Notation -
-
-

-
-
State
-

- A type that is the state_type - of the ODE -

-
WrappedState
-

- A type that is a model of State Wrapper for the state type State. -

-
x
-

- Object of type State -

-
w
-

- Object of type WrappedState -

-
-
-
- - Valid - Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Get resizeability -

- 		-

- is_resizeable< - State > -

- 		-

- boost::false_type or boost::true_type -

- 		-

- Returns boost::true_type if the State is resizeable, boost::false_type otherwise. -

-
-

- Create WrappedState - type -

- 		-

- state_wrapper< - State > -

- 		-

- WrappedState -

- 		-

- Creates the type for a WrappedState - for the state type State -

-
-

- Constructor -

- 		-

- WrappedState() -

- 		-

- WrappedState -

- 		-

- Constructs a state wrapper with an empty state -

-
-

- Copy Constructor -

- 		-

- WrappedState( - w ) -

- 		-

- WrappedState -

- 		-

- Constructs a state wrapper with a state of the same size as the - state in w -

-
-

- Get state -

- 		-

- w.m_v -

- 		-

- State -

- 		-

- Returns the State - object of this state wrapper. -

-
-

- Check size -

- 		-

- w.same_size( - x ) -

- 		-

- bool -

- 		-

- Returns true if x and w.m_v - have the same size, false - otherwise. Only required if State is resizeable. -

-
-

- Resize -

- 		-

- w.resize( - x ) -

- 		-

- bool -

- 		-

- If w.m_v and x - have different sizes, w.m_v - is resized to the size of x - and the method returns true. - Otherwise the method returns false. - Only required if State - is resizeable. -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHome -
- - diff --git a/doc/boost_numeric_odeint/concepts/stepper.html b/doc/boost_numeric_odeint/concepts/stepper.html deleted file mode 100644 index 43ae063c..00000000 --- a/doc/boost_numeric_odeint/concepts/stepper.html +++ /dev/null @@ -1,271 +0,0 @@ - - - -Stepper - - - - - - - - -
-PrevUpHomeNext -
-
-

-Stepper -

-

- This concepts specifies the interface a simple stepper has to fulfill to - be used within the integrate - functions. -

-
- - Description -
-

- The basic stepper concept. A basic stepper following this Stepper concept - is able to perform a single step of the solution x(t) - of an ODE to obtain x(t+dt) using a given step size - dt. Basic steppers can be Runge Kutta steppers, symplectic - steppers as well as implicit steppers. Depending on the actual stepper, the - ODE is defined as System, - Symplectic - System, Simple - Symplectic System or Implicit - System. Note that all error steppers are also basic steppers. -

-
- - Refinement - of -
-
-
- - Associated - types -
-
-
- - Notation -
-
-

-
-
Stepper
-

- A type that is a model of Stepper -

-
State
-

- A type representing the state x of the ODE -

-
Time
-

- A type representing the time t of the ODE -

-
stepper
-

- An object of type Stepper -

-
x
-

- Object of type State -

-
t, dt
-

- Objects of type Time -

-
sys
-

- An object defining the ODE. Depending on the Stepper this might be - a model of System, - Symplectic - System, Simple - Symplectic System or Implicit - System -

-
-
-
- - Valid Expressions -
-
------ - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Get the order -

- 		-

- stepper.order() -

- 		-

- order_type -

- 		-

- Returns the order of the stepper. -

-
-

- Do step -

- 		-

- stepper.do_step( - sys , - x , - t , - dt ) -

- 		-

- void -

- 		-

- Performs one step of step size dt. - The newly obtained state is written in place in x. -

-
-
- - Models -
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/symplectic_system.html b/doc/boost_numeric_odeint/concepts/symplectic_system.html deleted file mode 100644 index 922a7e61..00000000 --- a/doc/boost_numeric_odeint/concepts/symplectic_system.html +++ /dev/null @@ -1,233 +0,0 @@ - - - -Symplectic System - - - - - - - - -
-PrevUpHomeNext -
-
-

-Symplectic - System -

-
- - Description -
-

- This concept describes how to define a symplectic system written with generalized - coordinate q and generalized - momentum p: -

-

- q'(t) = f(p) -

-

- p'(t) = g(q) -

-

- Such a situation is typically found for Hamiltonian systems with a separable - Hamiltonian: -

-

- H(p,q) = Hkin(p) + V(q) -

-

- which gives the equations of motion: -

-

- q'(t) = dHkin / dp = f(p) -

-

- p'(t) = dV / dq = g(q) -

-

- The algorithmic implementation of this situation is described by a pair of - callable objects for f and g with - a specific parameter signature. Such a system should be implemented as a - std::pair of functions or a functors. Symplectic systems are used in symplectic - steppers like symplectic_rkn_sb3a_mclachlan. -

-
- - Notation -
-
-

-
-
System
-

- A type that is a model of Symplectic_System -

-
Coor
-

- The type of the coordinate q -

-
Momentum
-

- The type of the momentum p -

-
CoorDeriv
-

- The type of the derivative of coordinate q' -

-
MomentumDeriv
-

- The type of the derivative of momentum p' -

-
sys
-

- An object of the type System -

-
q
-

- Object of type Coor -

-
p
-

- Object of type Momentum -

-
dqdt
-

- Object of type CoorDeriv -

-
dpdt
-

- Object of type MomentumDeriv -

-
-
-
- - Valid - expressions -
-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Check for pair -

- 		-

- boost::is_pair< - System >::type -

- 		-

- boost::mpl::true_ -

- 		-

- Check if System is a pair -

-
-

- Calculate dq/dt = f(p) -

- 		-

- sys.first( - p , - dqdt ) -

- 		-

- void -

- 		-

- Calculates f(p), the result is stored into - dqdt -

-
-

- Calculate dp/dt = g(q) -

- 		-

- sys.second( - q , - dpdt ) -

- 		-

- void -

- 		-

- Calculates g(q), the result is stored into - dpdt -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/concepts/system.html b/doc/boost_numeric_odeint/concepts/system.html deleted file mode 100644 index ccbf96e7..00000000 --- a/doc/boost_numeric_odeint/concepts/system.html +++ /dev/null @@ -1,150 +0,0 @@ - - - -System - - - - - - - - -
-PrevUpHomeNext -
-
-

-System -

-
- - Description -
-

- The System concept models the algorithmic implementation of the rhs. of the - ODE x' = f(x,t). The only requirement for this concept - is that it should be callable with a specific parameter syntax (see below). - A System is typically implemented as a function or a functor. Systems fulfilling - this concept are required by all Runge-Kutta steppers as well as the Bulirsch-Stoer - steppers. However, symplectic and implicit steppers work with other system - concepts, see Symplectic - System and Implicit - System. -

-
- - Notation -
-
-

-
-
System
-

- A type that is a model of System -

-
State
-

- A type representing the state x of the ODE -

-
Deriv
-

- A type representing the derivative x' of the ODE -

-
Time
-

- A type representing the time -

-
sys
-

- An object of type System -

-
x
-

- Object of type State -

-
dxdt
-

- Object of type Deriv -

-
t
-

- Object of type Time -

-
-
-
- - Valid expressions -
-
------ - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Calculate dx/dt := f(x,t) -

- 		-

- sys( - x , - dxdt , - t ) -

- 		-

- void -

- 		-

- Calculates f(x,t), the result is stored into dxdt -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/getting_started.html b/doc/boost_numeric_odeint/getting_started.html deleted file mode 100644 index fd0d55e4..00000000 --- a/doc/boost_numeric_odeint/getting_started.html +++ /dev/null @@ -1,42 +0,0 @@ - - - -Getting started - - - - - - - - -
-PrevUpHomeNext -
-
-

-Getting started -

-
-
Overview
-
Usage, - Compilation, Headers
-
Short - Example
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/getting_started/overview.html b/doc/boost_numeric_odeint/getting_started/overview.html deleted file mode 100644 index 0666bdeb..00000000 --- a/doc/boost_numeric_odeint/getting_started/overview.html +++ /dev/null @@ -1,1114 +0,0 @@ - - - -Overview - - - - - - - - -
-PrevUpHomeNext -
-
-

-Overview -

-
- - - - - -
[Caution]Caution

- Boost.Numeric.Odeint is not an official boost library! -

-

- odeint is a library for solving initial value problems (IVP) of ordinary - differential equations. Mathematically, these problems are formulated as - follows: -

-

- x'(t) = f(x,t), x(0) = x0. -

-

- 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. -

-

- Ordinary 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 - in one coordinate. Then, a system of coupled ordinary differential occurs, - sometimes also referred as lattices ODEs. -

-

- Numerical 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 ~ dt2. 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. -

-

- The 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. In 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: -

-
-

- In odeint, the following algorithms are implemented: -

-
-

Table 1.1. Stepper Algorithms

-
----------- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Algorithm -

- 		-

- Class -

- 		-

- Concept -

- 		-

- System Concept -

- 		-

- Order -

- 		-

- Error Estimation -

- 		-

- Dense Output -

- 		-

- Internal state -

- 		-

- Remarks -

-
-

- Explicit Euler -

- 		-

- euler -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- 1 -

- 		-

- No -

- 		-

- Yes -

- 		-

- No -

- 		-

- Very simple, only for demonstrating purpose -

-
-

- Modified Midpoint -

- 		-

- modified_midpoint -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable (2) -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Used in Bulirsch-Stoer implementation -

-
-

- Runge-Kutta 4 -

- 		-

- runge_kutta4 -

- 		-

- Stepper -

- 		-

- System -

- 		-

- 4 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- The classical Runge Kutta scheme, good general scheme without error - control -

-
-

- Cash-Karp -

- 		-

- runge_kutta_cash_karp54 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 5 -

- 		-

- Yes (4) -

- 		-

- No -

- 		-

- No -

- 		-

- Good general scheme with error estimation, to be used in controlled_error_stepper -

-
-

- Dormand-Prince 5 -

- 		-

- runge_kutta_dopri5 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 5 -

- 		-

- Yes (4) -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- Standard method with error control and dense output, to be used - in controlled_error_stepper and in dense_output_controlled_explicit_fsal. -

-
-

- Fehlberg 78 -

- 		-

- runge_kutta_fehlberg78 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 8 -

- 		-

- Yes (7) -

- 		-

- No -

- 		-

- No -

- 		-

- Good high order method with error estimation, to be used in controlled_error_stepper. -

-
-

- Adams Bashforth -

- 		-

- adams_bashforth -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Multistep method -

-
-

- Adams Moulton -

- 		-

- adams_moulton -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Multistep method -

-
-

- Adams Bashforth Moulton -

- 		-

- adams_bashforth_moulton -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Combined multistep method -

-
-

- Controlled Runge Kutta -

- 		-

- controlled_runge_kutta -

- 		-

- Controlled - Stepper -

- 		-

- System -

- 		-

- depends -

- 		-

- Yes -

- 		-

- No -

- 		-

- depends -

- 		-

- Error control for Error - Stepper. Requires an Error - Stepper from above. Order depends on the given ErrorStepper -

-
-

- Dense Output Runge Kutta -

- 		-

- dense_output_runge_kutta -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- depends -

- 		-

- No -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- Dense output for Stepper - and Error - Stepper from above if they provide dense output functionality - (like euler and - runge_kutta_dopri5). - Order depends on the given stepper. -

-
-

- Bulirsch-Stoer -

- 		-

- bulirsch_stoer -

- 		-

- Controlled - Stepper -

- 		-

- System -

- 		-

- variable -

- 		-

- Yes -

- 		-

- No -

- 		-

- No -

- 		-

- Stepper with step size and order control. Very good if high precision - is required. -

-
-

- Bulirsch-Stoer Dense Output -

- 		-

- bulirsch_stoer_dense_out -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- variable -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Stepper with step size and order control as well as dense output. - Very good if high precision and dense output is required. -

-
-

- Implicit Euler -

- 		-

- implicit_euler -

- 		-

- Stepper -

- 		-

- Implicit - System -

- 		-

- 1 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Basic implicit routine. Requires the Jacobian. Works only with - Boost.UBlas - vectors as state types. -

-
-

- Rosenbrock 4 -

- 		-

- rosenbrock4 -

- 		-

- Error - Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Good for stiff systems. Works only with Boost.UBlas - vectors as state types. -

-
-

- Controlled Rosenbrock 4 -

- 		-

- rosenbrock4_controller -

- 		-

- Controlled - Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Rosenbrock 4 with error control. Works only with Boost.UBlas - vectors as state types. -

-
-

- Dense Output Rosenbrock 4 -

- 		-

- rosenbrock4_dense_output -

- 		-

- Dense - Output Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Controlled Rosenbrock 4 with dense output. Works only with Boost.UBlas - vectors as state types. -

-
-

- Symplectic Euler -

- 		-

- symplectic_euler -

- 		-

- Stepper -

- 		-

- Symplectic - System Simple - Symplectic System -

- 		-

- 1 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Basic symplectic solver for separable Hamiltonian system -

-
-

- Symplectic RKN McLachlan -

- 		-

- symplectic_rkn_sb3a_mclachlan -

- 		-

- Stepper -

- 		-

- Symplectic - System Simple - Symplectic System -

- 		-

- 6 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Symplectic solver for separable Hamiltonian system with order 6 -

-
-
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/getting_started/short_example.html b/doc/boost_numeric_odeint/getting_started/short_example.html deleted file mode 100644 index 2f9e5a2d..00000000 --- a/doc/boost_numeric_odeint/getting_started/short_example.html +++ /dev/null @@ -1,200 +0,0 @@ - - - -Short Example - - - - - - - - -
-PrevUpHomeNext -
-
-

-Short - Example -

-

- Imaging, you want to numerically integrate a harmonic oscillator with friction. - The equations of motion are given by x'' = -x + γ x'. - Odeint only deals with first order ODEs that have no higher derivatives than - x' involved. However, any higher order ODE can be transformed to a system - of first 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+γ p): -

-

-

-
/* The type of container used to hold the state vector */
-typedef std::vector< double > state_type;
-
-const double gam = 0.15;
-
-/* The rhs of x' = f(x) */
-void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
-{
-    dxdt[0] = x[1];
-    dxdt[1] = -x[0] - gam*x[1];
-}
-
-

-

-

- Here we chose vector<double> - as the state type, but others are also possible, for example boost::array<double,2>. odeint is designed in such a way that - you can easily use your own state types. Next, the ODE is defined which is - in this case a simple function calculating f(x)'. The - parameter signature of this function is crucial: the integration methods - will always call 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. -

-

- Now, we have to define the initial state from which the integration should - start: -

-

-

-
state_type x(2);
-x[0] = 1.0; // start at x=1.0, p=0.0
-x[1] = 0.0;
-
-

-

-

- For 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_kutta_rk5_ck - stepper (5th order) and uses adaptive step-size. -

-

-

-
size_t steps = integrate( harmonic_oscillator ,
-        x , 0.0 , 10.0 , 0.1 );
-
-

-

-

- The integrate function expects as parameters the rhs of the ode as defined - above, the initial state x, - the start-and end-time of the integration as well as the initial time step=size. - Note, that integrate - uses an adaptive step-size during the integration steps so the time points - will not be equally spaced. The integration returns the number of steps that - were applied and updates x which is set to the approximate solution of the - ODE at the end of integration. -

-

- It is, of course, also possible to represent the ode system as a class. The - rhs must then be implemented as a functor having defined the ()-operator: -

-

-

-
/* The rhs of x' = f(x) defined as a class */
-class harm_osc {
-
-    double m_gam;
-
-public:
-    harm_osc( double gam ) : m_gam(gam) { }
-
-    void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
-    {
-        dxdt[0] = x[1];
-        dxdt[1] = -x[0] - m_gam*x[1];
-    }
-};
-
-

-

-

- which can be used via -

-

-

-
harm_osc ho(0.15);
-steps = integrate( ho ,
-        x , 0.0 , 10.0 , 0.1 );
-
-

-

-

- You surely have already noticed that during the integration a lot of steps - had to be done. You might wonder if you could access them do observe the - solution during the iteration. Yes, this is possible, of course. All you - have to do is to provide a reasonable observer. An example is -

-

-

-
struct push_back_state_and_time
-{
-    std::vector< state_type >& m_states;
-    std::vector< double >& m_times;
-
-    push_back_state_and_time( std::vector< state_type > &states , std::vector< double > &times )
-    : m_states( states ) , m_times( times ) { }
-
-    void operator()( const state_type &x , double t )
-    {
-        m_states.push_back( x );
-        m_times.push_back( t );
-    }
-};
-
-

-

-

- which stores the intermediate steps in a container. Note, the argument structure - of the ()-operator: odeint calls the observer exactly in this way, providing - the current state and time. Now, you only have to pass this container to - the integration function: -

-

-

-
vector<state_type> x_vec;
-vector<double> times;
-
-steps = integrate( harmonic_oscillator ,
-        x , 0.0 , 10.0 , 0.1 ,
-        push_back_state_and_time( x_vec , times ) );
-
-/* output */
-for( size_t i=0; i<=steps; i++ )
-{
-    cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
-}
-
-

-

-

- That is all. Of course, you can use functional libraries like Boost.Lambda - or Boost.Phoenix - to ease the creation of observer functions. -

-

- The full cpp file for this example can be found here: ../../examples/harmonic_oscillator.cpp -

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/getting_started/usage__compilation__headers.html b/doc/boost_numeric_odeint/getting_started/usage__compilation__headers.html deleted file mode 100644 index 8275fd09..00000000 --- a/doc/boost_numeric_odeint/getting_started/usage__compilation__headers.html +++ /dev/null @@ -1,52 +0,0 @@ - - - -Usage, Compilation, Headers - - - - - - - - -
-PrevUpHomeNext -
-
-

-Usage, - Compilation, Headers -

-

- odeint is completely header-only, no linking against pre-compiled code is - required. It can be include by -

-

-

-
#include <boost/numeric/odeint.hpp>
-
-

- which includes all headers of the library. All functions and classes from - odeint live in the namespace -

-
using namespace boost::numeric::odeint;
-
-

-

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/odeint_in_detail.html b/doc/boost_numeric_odeint/odeint_in_detail.html deleted file mode 100644 index 8a40dea7..00000000 --- a/doc/boost_numeric_odeint/odeint_in_detail.html +++ /dev/null @@ -1,49 +0,0 @@ - - - -odeint in detail - - - - - - - - -
-PrevUpHomeNext -
-
-

-odeint in detail -

-
-
Steppers
-
Generation - functions
-
Integrate - functions
-
Iterators
-
State - types, algebras and operations
-
Using - boost::ref
-
Using - boost::range
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/odeint_in_detail/generation_functions.html b/doc/boost_numeric_odeint/odeint_in_detail/generation_functions.html deleted file mode 100644 index 573fd9ea..00000000 --- a/doc/boost_numeric_odeint/odeint_in_detail/generation_functions.html +++ /dev/null @@ -1,301 +0,0 @@ - - - -Generation functions - - - - - - - - -
-PrevUpHomeNext -
-
-

-Generation - functions -

-

- In the Tutorial we have - learned how we can use the generation functions make_controlled - and make_dense_output to - create controlled and dense output stepper from a simple stepper or an error - stepper. The syntax of these two functions is very simple: -

-

-

-
auto stepper1 = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() );
-auto stepper2 = make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() );
-
-

-

-

- The first two parameters are the absolute and the relative error tolerances - and the third parameter is the stepper. In C++03 you can infer the type from - the result_of mechanism: -

-

-

-
boost::numeric::odeint::result_of::make_controlled< stepper_type >::type stepper3 = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() );
-boost::numeric::odeint::result_of::make_dense_output< stepper_type >::type stepper4 = make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() );
-
-

-

-

- To use your own steppers with the make_controlled - or make_dense_output you - need to specialize two class templates. Suppose your steppers are called - custom_stepper, custom_controller and custom_dense_output. - Then, the first class you need to specialize is boost::numeric::get_controller, - a meta function returning the type of the controller: -

-

-

-
namespace boost { namespace numeric { namespace odeint {
-
-template<>
-struct get_controller< custom_stepper >
-{
-    typedef custom_controller type;
-};
-
-} } }
-
-

-

-

- The second one is a factory class boost::numeric::odeint::controller_factory - which constructs the controller from the tolerances and the stepper. In our - dummy implementation this class is -

-

-

-
namespace boost { namespace numeric { namespace odeint {
-
-template<>
-struct controller_factory< custom_stepper , custom_controller >
-{
-    custom_controller operator()( double abs_tol , double rel_tol , const custom_stepper & ) const
-    {
-        return custom_controller();
-    }
-};
-
-} } }
-
-

-

-

- This is all to use the make_controlled - mechanism. Now you can use your controller via -

-

-

-
auto stepper5 = make_controlled( 1.0e-6 , 1.0e-6 , custom_stepper() );
-
-

-

-

- For the dense_output_stepper everything works similar. Here you have to specialize - boost::numeric::odeint::get_dense_output and boost::numeric::odeint::dense_output_factory. - These two classes have the same syntax as their relatives get_controller - and controller_factory. -

-

- Of course, all controllers and dense-output steppers in odeint can be used - with these mechanisms. In the table below you will find, which steppers is - constructed from make_controlled - or make_dense_output if applied - on a stepper from odeint: -

-
-

Table 1.8. Generation functions make_controlled( abs_error , rel_error , stepper - )

-
----- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Stepper -

- 		-

- Result of make_controlled -

- 		-

- Remarks -

-
-

- runge_kutta_cash_karp54 -

- 		-

- controlled_runge_kutta< runge_kutta_cash_karp54 - , default_error_checker<...> > -

- 		-

- ax=1, adxdt=1 -

-
-

- runge_kutta_fehlberg78 -

- 		-

- controlled_runge_kutta< runge_kutta_fehlberg78 - , default_error_checker<...> > -

- 		-

- ax=1, adxdt=1 -

-
-

- runge_kutta_dopri5 -

- 		-

- controlled_runge_kutta< runge_kutta_dopri5 - , default_error_checker<...> > -

- 		-

- a x=1, adxdt=1 -

-
-

- rosenbrock4 -

- 		-

- rosenbrock4_controlled< rosenbrock4 - > -

- 		-

- - -

-
-
-
-

Table 1.9. Generation functions make_dense_output( abs_error , rel_error , stepper - )

-
----- - - - - - - - - - - - - - - - - - -
-

- Stepper -

- 		-

- Result of make_dense_output -

- 		-

- Remarks -

-
-

- runge_kutta_dopri5 -

- 		-

- dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5 - , default_error_checker<...> > - > -

- 		-

- a x=1, adxdt=1 -

-
-

- rosenbrock4 -

- 		-

- rosenbrock4_dense_output< rosenbrock4_controller< rosenbrock4 - > > -

- 		-

- - -

-
-
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/odeint_in_detail/integrate_functions.html b/doc/boost_numeric_odeint/odeint_in_detail/integrate_functions.html deleted file mode 100644 index 057c6522..00000000 --- a/doc/boost_numeric_odeint/odeint_in_detail/integrate_functions.html +++ /dev/null @@ -1,340 +0,0 @@ - - - -Integrate functions - - - - - - - - -
-PrevUpHomeNext -
-
-

-Integrate - functions -

-

- Integrate functions perform the time evolution of a given ODE from some starting - time t0 to a given end time t1 - and starting at state x0 by subsequent calls of a given - stepper's do_step function. - Additionally, the user can provide an __observer to analyze the state during - time evolution. There are five different integrate functions which have different - strategies on when to call the observer function during integration. All - of the integrate functions except integrate_n_steps - can be called with any stepper following one of the stepper concepts: Stepper , Error - Stepper , Controlled - Stepper , Dense - Output Stepper. Depending on the abilities of the stepper, the integrate - functions make use of step-size control or dense output. -

-
- - Equidistant - observer calls -
-

- If observer calls at equidistant time intervals dt are - needed, the integrate_const - or integrate_n_steps function - should be used. We start with explaining integrate_const: -

-

- integrate_const( - stepper , - system , - x0 , - t0 , - t1 , - dt ) -

-

- integrate_const( - stepper , - system , - x0 , - t0 , - t1 , - dt , - observer ) -

-

- These integrate the ODE given by system - with subsequent steps from stepper. - Integration start at t0 and - x0 and ends at some t' - = t0 + n dt with n such that t1 - - dt < t' <= t1. x0 - is changed to the approximative solution x(t') at the - end of integration. If provided, the observer - is invoked at times t0, t0 + dt, - t0 + 2dt, ... ,t'. integrate_const returns the number of steps - performed during the integration. Note that if you are using a simple Stepper or Error - Stepper and want to make exactly n - steps you should prefer the integrate_n_steps - function below. -

-
-
- - Integrate - a given number of steps -
-

- This function is very similar to integrate_const - above. The only difference is that it does not take the end time as parameter, - but rather the number of steps. The integration is then performed until the - time t0+n*dt. -

-

- integrate_n_steps( - stepper , - system , - x0 , - t0 , - dt , - n ) -

-

- integrate_n_steps( - stepper , - system , - x0 , - t0 , - dt , - n , observer ) -

-

- Integrates the ODE given by system - with subsequent steps from stepper - starting at x0 and t0. If provided, - observer is called after - every step and at the beginning with t0, - similar as above. The approximate result for x( t0 + n dt ) - is stored in x0. This function - returns the end time t0 + n*dt. -

-
- - Observer - calls at each step -
-

- If the observer should be called at each time step then the integrate_adaptive function should be used. - Note that in the case of Controlled - Stepper or Dense - Output Stepper this leads to non-equidistant observer calls as the - step size changes. -

-

- integrate_adaptive( - stepper , - system , - x0 , - t0 , - t1 , - dt ) -

-

- integrate_adaptive( - stepper , - system , - x0 , - t0 , - t1 , - dt , - observer ) -

-

- Integrates the ODE given by system - with subsequent steps from stepper. - Integration start at t0 and - x0 and ends at t1. - x0 is changed to the approximative - solution x(t1) at the end of integration. If provided, - the observer is called after - each step (and before the first step at t0). - integrate_adaptive returns - the number of steps performed during the integration. -

-
-
- - Observer - calls at given time points -
-

- If the observer should be called at some user given time points the integrate_times function should be used. - The times for observer calls are provided as a sequence of time values. The - sequence is either defined via two iterators pointing to begin and end of - the sequence or in terms of a Boost.Range - object. -

-

- integrate_times( - stepper , - system , - x0 , - times_start , - times_end , - dt , - observer ) -

-

- integrate_times( - stepper , - system , - x0 , - time_range , - dt , - observer ) -

-

- Integrates the ODE given by system - with subsequent steps from stepper. - Integration starts at *times_start - and ends exactly at *(times_end-1). - x0 contains the approximate - solution at the end point of integration. This function requires an observer - which is invoked at the subsequent times *times_start++ - until times_start == - times_end. If called with a Boost.Range - time_range the function behaves - the same with times_start = boost::begin( time_range - ) and times_end - = boost::end( - time_range ). - integrate_times returns the - number of steps performed during the integration. -

-
-
- - Convenience - integrate function -
-

- Additionally to the sophisticated integrate function above odeint also provides - a simple integrate routine - which uses a dense output stepper based on runge_kutta_dopri5 - with standard error bounds 10-6 for the steps. -

-

- integrate( - system , - x0 , - t0 , - t1 , - dt ) -

-

- integrate( - system , - x0 , - t0 , - t1 , - dt , - observer ) -

-

- This function behaves exactly like integrate_adaptive - above but no stepper has to be provided. It also returns the number of steps - performed during the integration. -

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html b/doc/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html deleted file mode 100644 index aae3772b..00000000 --- a/doc/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html +++ /dev/null @@ -1,1327 +0,0 @@ - - - -State types, algebras and operations - - - - - - - - -
-PrevUpHomeNext -
-
-

-State - types, algebras and operations -

-
-
Construction/Resizing
-
Algebras - and Operations
-
Adapt - your own operations
-
-

- In odeint the stepper algorithms are implemented independently of the underlying - fundamental mathematical operations. This is realized by giving the user - full control over the state type and the mathematical operations for this - state type. Technically, this is done by introducing three concepts: StateType, - Algebra, Operations. Most of the steppers in odeint expect three class types - fulfilling these concepts as template parameters. Note that these concepts - are not fully independent of each other but rather a valid combination must - be provided in order to make the steppers work. In the following we will - give some examples on reasonable state_type-algebra-operations combinations. - For the most common state types, like vector<double> or array<double,N> - the default values range_algebra and default_operations are perfectly fine - and odeint can be used as is without worrying about algebra/operations at - all. -

-
- - - - - -
[Important]Important

- state_type, algebra and operations are not independent, a valid combination - must be provided to make odeint work properly -

-

- Moreover, as odeint handles the memory required for intermediate temporary - objects itself, it also needs knowledge about how to create state_type objects - and maybe how to allocate memory (resizing). All in all, the following things - have to be taken care of when odeint is used with non-standard state types: -

-
-

- Again, odeint already provides basic interfaces for most of the usual state - types. So if you use a std::vector, - or a boost::array as state type no additional work - is required, they just work out of the box. -

-
-

-Construction/Resizing -

-
-
Using - the container interface
-
std::list
-
-

- We distinguish between two basic state types: fixed sized and dynamically - sized. For fixed size state types the default constructor state_type() - already allocates the required memory, prominent example is boost::array<T,N>. Dynamically sized types have to be - resized to make sure enough memory is allocated, the standard constructor - does not take care of the resizing. Examples for this are the STL containers - like vector<double>. -

-

- The most easy way of getting your own state type to work with odeint is - to use a fixed size state, base calculations on the range_algebra and provide - the following functionality: -

-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Construct State -

- 		-

- State x() -

- 		-

- void -

- 		-

- Creates an instance of State - and allocates memory. -

-
-

- Begin of the sequence -

- 		-

- boost::begin(x) -

- 		-

- Iterator -

- 		-

- Returns an iterator pointing to the begin of the sequence -

-
-

- End of the sequence -

- 		-

- boost::end(x) -

- 		-

- Iterator -

- 		-

- Returns an iterator pointing to the end of the sequence -

-
-
- - - - - -
[Warning]Warning

- If your state type does not allocate memory by default construction, - you must define it as resizeable and - provide resize functionality (see below). Otherwise segmentation faults - will occur. -

-

- So fixed sized arrays supported by Boost.Range - immediately work with odeint. For dynamically sized arrays one has to additionally - supply the resize functionality. First, the state has to be tagged as resizeable - by specializing the struct is_resizeable - which consists of one typedef and one bool value: -

-
------ - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Resizability -

- 		-

- is_resizeable<State>::type -

- 		-

- boost::true_type or boost::false_type -

- 		-

- Determines resizeability of the state type, returns boost::true_type if the state is resizeable. -

-
-

- Resizability -

- 		-

- is_resizeable<State>::value -

- 		-

- bool -

- 		-

- Same as above, but with bool - value. -

-
-

- Defining type to be true_type and value - as true tells odeint that - your state is resizeable. By default, odeint now expects the support of - boost::size(x) and - a x.resize( boost::size(y) ) - member function for resizing: -

-
------ - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Get size -

- 		-

- boost::size( - x ) -

- 		-

- size_type -

- 		-

- Returns the current size of x. -

-
-

- Resize -

- 		-

- x.resize( - boost::size( - y ) - ) -

- 		-

- void -

- 		-

- Resizes x to have the same size as y. -

-
-
-
-

- As a first example we take the most simple case and implement our own - vector my_vector which - will provide a container interface. This makes Boost.Range - working out-of-box. We add a little functionality to our vector which - makes it allocate some default capacity by construction. This is helpful - when using resizing as then a resize can be assured to not require a - new allocation. -

-

-

-
template< int MAX_N >
-class my_vector
-{
-    typedef std::vector< double > vector;
-
-public:
-    typedef vector::iterator iterator;
-    typedef vector::const_iterator const_iterator;
-
-public:
-    my_vector( const size_t N )
-        : m_v( N )
-    {
-        m_v.reserve( MAX_N );
-    }
-
-    my_vector()
-        : m_v()
-    {
-        m_v.reserve( MAX_N );
-    }
-
-// ... [ implement container interface ]
-
-

-

-

- The only thing that has to be done other than defining is thus declaring - my_vector as resizeable: -

-

-

-
// define my_vector as reizeable
-
-namespace boost { namespace numeric { namespace odeint {
-
-template<size_t N>
-struct is_resizeable< my_vector<N> >
-{
-    typedef boost::true_type type;
-    static const bool value = type::value;
-};
-
-} } }
-
-

-

-

- If we wouldn't specialize the is_resizeable - template, the code would still compile but odeint would not adjust the - size of temporary internal instances of my_vector and hence try to fill - zero-sized vectors resulting in segmentation faults! The full example - can be found in my_vector.cpp -

-
-
-
-

- If your state type does work with Boost.Range, - but handles resizing differently you are required to specialize two implementations - used by odeint to check a state's size and to resize: -

-
------ - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Check size -

- 		-

- same_size_impl<State,State>::same_size(x - , y) -

- 		-

- bool -

- 		-

- Returns true if the size of x equals the size of y. -

-
-

- Resize -

- 		-

- resize_impl<State,State>::resize(x - , y) -

- 		-

- void -

- 		-

- Resizes x to have the same size as y. -

-
-

- As an example we will use a std::list - as state type in odeint. Because std::list - is not supported by boost::size - we have to replace the same_size and resize implementation to get list - to work with odeint. The following code shows the required template specializations: -

-

-

-
typedef std::list< double > state_type;
-
-namespace boost { namespace numeric { namespace odeint {
-
-template< >
-struct is_resizeable< state_type >
-{ // declare resizeablility
-    typedef boost::true_type type;
-    const static bool value = type::value;
-};
-
-template< >
-struct same_size_impl< state_type , state_type >
-{ // define how to check size
-    static bool same_size( const state_type &v1 ,
-                           const state_type &v2 )
-    {
-        return v1.size() == v2.size();
-    }
-};
-
-template< >
-struct resize_impl< state_type , state_type >
-{ // define how to resize
-    static void resize( state_type &v1 ,
-                        const state_type &v2 )
-    {
-        v1.resize( v2.size() );
-    }
-};
-
-} } }
-
-

-

-

- With these definitions odeint knows how to resize std::lists - and so they can be used as state types. A complete example can be found - in list_lattice.cpp. -

-
-
-
-

-Algebras - and Operations -

-
-
GSL - Vector
-
Vector - Space Algebra
-
Boost.Ublas
-
Point - type
-
-

- To provide maximum flexibility odeint is implemented in a highly modularized - way. This means it is possible to change the underlying mathematical operations - without touching the integration algorithms. The fundamental mathematical - operations are those of a vector space, that is addition of state_types and multiplication of state_types with a scalar (time_type). In odeint this is realized - in two concepts: Algebra and Operations. The standard way how this works - is by the range algebra which provides functions that apply a specific - operation to each of the individual elements of a container based on the - Boost.Range - library. If your state type is not supported by Boost.Range - there are several possibilities to tell odeint how to do algebraic operations: -

-
    -
  • - Implement boost::begin and boost::end - for your state type so it works with Boost.Range. -
    • -
  • - Implement vector-vector addition operator + - and scalar-vector multiplication operator * - and use the non-standard vector_space_algebra. -
    • -
  • - Implement your own algebra that implements the required functions. -
    • -
-
-
-

- In the following example we will try to use the gsl_vector - type from GSL (GNU Scientific - Library) as state type in odeint. We will realize this by implementing - a wrapper around the gsl_vector that takes care of construction/destruction. - Also, Boost.Range - is extended such that it works with gsl_vectors - as well which required also the implementation of a new gsl_iterator. -

-
- - - - - -
[Note]Note

- odeint already includes all the code presented here, see gsl_wrapper.hpp, - so gsl_vectors can - be used straight out-of-box. The following description is just for - educational purpose. -

-

- The GSL is a C library, so gsl_vector - has neither constructor, nor destructor or any begin - or end function, no iterators - at all. So to make it work with odeint plenty of things have to be implemented. - Note that all of the work shown here is already included in odeint, so - using gsl_vectors in - odeint doesn't require any further adjustments. We present it here just - as an educational example. We start with defining appropriate constructors - and destructors. This is done by specializing the state_wrapper - for gsl_vector. State - wrappers are used by the steppers internally to create and manage temporary - instances of state types: -

-

-

-
template<>
-struct state_wrapper< gsl_vector* >
-{
-    typedef double value_type;
-    typedef gsl_vector* state_type;
-    typedef state_wrapper< gsl_vector* > state_wrapper_type;
-
-    state_type m_v;
-
-    state_wrapper( )
-    {
-        m_v = gsl_vector_alloc( 1 );
-    }
-
-    state_wrapper( const state_wrapper_type &x )
-    {
-        resize( m_v , x.m_v );
-        gsl_vector_memcpy( m_v , x.m_v );
-    }
-
-
-    ~state_wrapper()
-    {
-        gsl_vector_free( m_v );
-    }
-
-};
-
-

-

-

- This state_wrapper specialization - tells odeint how gsl_vectors are created, copied and destroyed. Next - we need resizing, this is required because gsl_vectors are dynamically - sized objects: -

-
template<>
-struct is_resizeable< gsl_vector* >
-{
-    typedef boost::true_type type;
-    const static bool value = type::value;
-};
-
-template <>
-struct same_size_impl< gsl_vector* , gsl_vector* >
-{
-    static bool same_size( const gsl_vector* x , const gsl_vector* y )
-    {
-        return x->size == y->size;
-    }
-};
-
-template <>
-struct resize_impl< gsl_vector* , gsl_vector* >
-{
-    static void resize( gsl_vector* x , const gsl_vector* y )
-    {
-        gsl_vector_free( x );
-        x = gsl_vector_alloc( y->size );
-    }
-};
-
-

-

-

- Up to now, we defined creation/destruction and resizing, but gsl_vectors - also don't support iterators, so we first implement a gsl iterator: -

-

-

-
/*
- * defines an iterator for gsl_vector
- */
-class gsl_vector_iterator
-      : public boost::iterator_facade< gsl_vector_iterator , double ,
-                                       boost::random_access_traversal_tag >
-{
-public :
-
-    gsl_vector_iterator( void ): m_p(0) , m_stride( 0 ) { }
-    explicit gsl_vector_iterator( gsl_vector *p ) : m_p( p->data ) , m_stride( p->stride ) { }
-    friend gsl_vector_iterator end_iterator( gsl_vector * );
-
-private :
-
-    friend class boost::iterator_core_access;
-    friend class const_gsl_vector_iterator;
-
-    void increment( void ) { m_p += m_stride; }
-    void decrement( void ) { m_p -= m_stride; }
-    void advance( ptrdiff_t n ) { m_p += n*m_stride; }
-    bool equal( const gsl_vector_iterator &other ) const { return this->m_p == other.m_p; }
-    bool equal( const const_gsl_vector_iterator &other ) const;
-    double& dereference( void ) const { return *m_p; }
-
-    double *m_p;
-    size_t m_stride;
-};
-
-

- A similar class exists for the const - version of the iterator. Then we have a function returning the end iterator - (similarly for const again): -

-
gsl_vector_iterator end_iterator( gsl_vector *x )
-{
-    gsl_vector_iterator iter( x );
-    iter.m_p += iter.m_stride * x->size;
-    return iter;
-}
-
-

-

-

- Finally, the bindings for Boost.Range - are added: -

-
// template<>
-inline gsl_vector_iterator range_begin( gsl_vector *x )
-{
-    return gsl_vector_iterator( x );
-}
-
-// template<>
-inline gsl_vector_iterator range_end( gsl_vector *x )
-{
-    return end_iterator( x );
-}
-
-

- Again with similar definitions for the const - versions. This eventually makes odeint work with gsl vectors as state - types. The full code for these bindings is found in gsl_wrapper.hpp. - It might look rather complicated but keep in mind that gsl is a pre-compiled - C library so it is quite an achievement that it's possible to get it - working at all in the first place. -

-
-
-
-

- As seen above, the standard way of performing algebraic operations on - container-like state types in odeint is to iterate through the elements - of the container and perform the operations element-wise on the underlying - value type. This is realized by means of the range_algebra - that uses Boost.Range - for obtaining iterators of the state types. However, there are other - ways to implement the algebraic operations on containers, one of which - is defining the addition/multiplication operators for the containers - directly and then using the vector_space_algebra. - If you use this algebra, the following operators have to be defined for - the state_type: -

-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Addition -

- 		-

- x + - y -

- 		-

- state_type -

- 		-

- Calculates the vector sum 'x+y'. -

-
-

- Assign addition -

- 		-

- x += - y -

- 		-

- state_type -

- 		-

- Performs x+y in place. -

-
-

- Scalar multiplication -

- 		-

- a * - x -

- 		-

- state_type -

- 		-

- Performs multiplication of vector x with scalar a. -

-
-

- Assign scalar multiplication -

- 		-

- x *= - a -

- 		-

- state_type -

- 		-

- Performs in-place multiplication of vector x with scalar a. -

-
-

- Defining these operators makes your state type work with any basic Runge-Kutta - stepper. However, if you want to use step-size control, some more functionality - is required. Specifically, operations like maxi( |erri| / (alpha - * |si|) ) have to be performed. err and - s are state_types, alpha is a scalar. As you can - see, we need element wise absolute value and division as well as an reduce - operation to get the maximum value. So for controlled steppers the following - things have to be implemented: -

-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Name -

- 		-

- Expression -

- 		-

- Type -

- 		-

- Semantics -

-
-

- Division -

- 		-

- x / - y -

- 		-

- state_type -

- 		-

- Calculates the element-wise division 'x/y' -

-
-

- Absolute value -

- 		-

- abs( - x ) -

- 		-

- state_type -

- 		-

- Element wise absolute value -

-
-

- Reduce -

- 		-

- vector_space_reduce_impl< state_type - >::reduce( state - , operator - , init - ) -

- 		-

- value_type -

- 		-

- Performs the operation operator - for subsequently each element of state - and returns the aggregate value. E.g. -

-

- init = - operator( - init , - state[0] - ); -

-

- init = - operator( - init , - state[1] - ) -

-

- ... -

-
-
-
-
-

- As an example for the employment of the vector_space_algebra - we will adopt ublas::vector from Boost.UBlas - to work as a state type in odeint. This is particularly easy because - ublas::vector supports vector-vector addition - and scalar-vector multiplication described above as well as boost::size. It also has a resize member function - so all that has to be done in this case is to declare resizability: -

-

-

-
typedef boost::numeric::ublas::vector< double > state_type;
-
-namespace boost { namespace numeric { namespace odeint {
-
-template<>
-struct is_resizeable< state_type >
-{
-    typedef boost::true_type type;
-    const static bool value = type::value;
-};
-
-} } }
-
-

-

-

- Now ublas::vector can be used as state type for simple Runge-Kutta steppers - in odeint by specifying the vector_space_algebra - as algebra in the template parameter list of the stepper. The following - code shows the corresponding definitions: -

-

-

-
int main()
-{
-    state_type x(3);
-    x[0] = 10.0; x[1] = 5.0 ; x[2] = 0.0;
-    typedef runge_kutta4< state_type , double , state_type , double , vector_space_algebra > stepper;
-    integrate_const( stepper() , lorenz , x ,
-                     0.0 , 10.0 , 0.1 );
-}
-
-

-

-

- Note again, that we haven't supported the requirements for controlled - steppers, but only for simple Runge-Kutta methods. You can find the full - example in lorenz_ublas.cpp. -

-
-
-
-

- Here we show how to implement the required operators on a state type. - As example we define a new class point3D - representing a three-dimensional vector with components x,y,z and define - addition and scalar multiplication operators for it. We use Boost.Operators - to reduce the amount of code to be written. The class for the point type - looks as follows: -

-

-

-
class point3D :
-    boost::additive1< point3D ,
-    boost::additive2< point3D , double ,
-    boost::multiplicative2< point3D , double > > >
-{
-public:
-
-    double x , y , z;
-
-    point3D()
-        : x( 0.0 ) , y( 0.0 ) , z( 0.0 )
-    { }
-
-    point3D( const double val )
-        : x( val ) , y( val ) , z( val )
-    { }
-
-    point3D( const double _x , const double _y , const double _z  )
-        : x( _x ) , y( _y ) , z( _z )
-    { }
-
-    point3D& operator+=( const point3D &p )
-    {
-        x += p.x; y += p.y; z += p.z;
-        return *this;
-    }
-
-    point3D& operator*=( const double a )
-    {
-        x *= a; y *= a; z *= a;
-        return *this;
-    }
-
-};
-
-

-

-

- By deriving from Boost.Operators - classes we don't have to define outer class operators like operator+( point3D , - point3D ) - because that is taken care of by the operators library. Note that for - simple Runge-Kutta schemes (like runge_kutta4) - only the + and * operators are required. If, however, - a controlled stepper is used one also needs to specify the division operator - / because calculation of - the error term involves an element wise division of the state types. - Additionally, controlled steppers require an abs - function calculating the element-wise absolute value for the state type: -

-

-

-
// only required for steppers with error control
-point3D operator/( const point3D &p1 , const point3D &p2 )
-{
-    return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p1.z );
-}
-
-point3D abs( const point3D &p )
-{
-    return point3D( std::abs(p.x) , std::abs(p.y) , std::abs(p.z) );
-}
-
-

-

-

- Finally, we have to add a specialization for reduce - implementing a reduction over the state type: -

-

-

-
namespace boost { namespace numeric { namespace odeint {
-// specialization of vector_space_reduce, only required for steppers with error control
-template<>
-struct vector_space_reduce< point3D >
-{
-    template< class Value , class Op >
-    Value operator() ( const point3D &p , Op op , Value init )
-    {
-        init = op( init , p.x );
-        //std::cout << init << " ";
-        init = op( init , p.y );
-        //std::cout << init << " ";
-        init = op( init , p.z );
-        //std::cout << init << std::endl;
-        return init;
-    }
-};
-} } }
-
-

-

-

- Again, note that the two last steps were only required if you want to - use controlled steppers. For simple steppers definition of the simple - += and *= - operators are sufficient. Having defined such a point type, we can easily - perform the integration on a Lorenz system by using the vector_space_algebra again: -

-

-

-
const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-void lorenz( const point3D &x , point3D &dxdt , const double t )
-{
-    dxdt.x = sigma * ( x.y - x.x );
-    dxdt.y = R * x.x - x.y - x.x * x.z;
-    dxdt.z = -b * x.z + x.x * x.y;
-}
-
-using namespace boost::numeric::odeint;
-
-int main()
-{
-
-    point3D x( 10.0 , 5.0 , 5.0 );
-    // point type defines it's own operators -> use vector_space_algebra !
-    typedef runge_kutta_dopri5< point3D , double , point3D ,
-                                double , vector_space_algebra > stepper;
-    int steps = integrate_adaptive( make_controlled<stepper>( 1E-10 , 1E-10 ) , lorenz , x ,
-                                    0.0 , 10.0 , 0.1 );
-    std::cout << x << std::endl;
-    std::cout << "steps: " << steps << std::endl;
-}
-
-

-

-

- The whole example can be found in lorenz_point.cpp -

-
-
-

- gsl_vector, gsl_matrix, ublas::matrix, blitz::matrix, thrust -

-
-

-Adapt - your own operations -

-

- to be continued -

-
    -
  • - thrust -
    • -
  • - gsl_complex -
    • -
  • - min, max, pow -
    • -
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/odeint_in_detail/steppers.html b/doc/boost_numeric_odeint/odeint_in_detail/steppers.html deleted file mode 100644 index 8b1bb6c8..00000000 --- a/doc/boost_numeric_odeint/odeint_in_detail/steppers.html +++ /dev/null @@ -1,2314 +0,0 @@ - - - -Steppers - - - - - - - - -
-PrevUpHomeNext -
-
-

-Steppers -

-
-
Explicit - steppers
-
Symplectic - solvers
-
Implicit - solvers
-
Multistep - methods
-
Controlled - steppers
-
Dense - output steppers
-
Using - steppers
-
Stepper - overview
-
Write - own steppers
-
Writing - own Runge-Kutta steppers
-
-

- Solving ordinary differential equation numerically is usually done iteratively, - that is a given state of an ordinary differential equation is iterated forward - x(t) -> x(t+dt) -> x(t+2dt). The steppers in odeint - perform one single step. The most general stepper type is described by the - Stepper concept. - The stepper concepts of odeint are described in detail in section Concepts, - here we briefly present the mathematical and numerical details of the steppers. - The Stepper - has two versions of the do_step - method, one with an in-place transform of the current state and one with - an out-of-place transform: -

-

- do_step( - sys , - inout , - t , dt ) -

-

- do_step( - sys , - in , - t , out , dt ) -

-

- The first parameter is always the system function - a function describing - the ODE. In the first version the second parameter is the step which is here - updated in-place and the third and the fourth parameters are the time and - step size (the time step). After a call of do_step - the state inout is updated - and now represents an approximate solution of the ODE at time t+dt. - In the second version the second argument is the state of the ODE at time - t, the third argument is t, the fourth argument is the - approximate solution at time t+dt which is filled by - do_step and the fifth argument - is the time step. Note that these functions do not change the time t. -

-

- System functions -

-

- Up to now, we have nothing said about the system function. This function - depends on the stepper. For the explicit Runge-Kutta steppers this function - can be a simple callable object hence a simple (global) C-function or a functor. - The parameter syntax is sys( x , - dxdt , - t ) - and it is assumed that it calculates dx/dt = f(x,t). - The function structure in most cases looks like: -

-

-

-
void sys( const state_type & /*x*/ , state_type & /*dxdt*/ , const double /*t*/ )
-{
-    // ...
-}
-
-

-

-

- Other types of system functions might represent Hamiltonian systems or system - which also compute the Jacobian needed in implicit steppers. For information - which stepper uses which system function see the stepper table below. It - might be possible that odeint will introduce new system types in near future. - Since the system function is strongly related to the stepper type, such an - introduction of a new stepper might result in a new type of system function. -

-
-

-Explicit - steppers -

-

- A first specialization are the explicit steppers. Explicit means that the - new state of the ode can be computed explicitly from the current state - without solving implicit equations. Such steppers have in common that they - evaluate the system at time t such that the result - of f(x,t) can be passed to the stepper. In odeint, - the explicit stepper have two additional methods -

-

- do_step( - sys , - inout , - dxdtin , - t , - dt ) -

-

- do_step( - sys , - in , - dxdtin , - t , - out , - dt ) -

-

- Here, the additional parameter is the value of the function f - at state x and time t. An example - is the Runge Kutta stepper of fourth order: -

-

-

-
runge_kutta4< state_type > rk;
-rk.do_step( sys1 , inout , t , dt );               // In-place transformation of inout
-rk.do_step( sys2 , inout , t , dt );               // call with different system: Ok
-rk.do_step( sys1 , in , t , out , dt );            // Out-of-place transformation
-rk.do_step( sys1 , inout , dxdtin , t , dt );      // In-place tranformation of inout
-rk.do_step( sys1 , in , dxdtin , t , out , dt );   // Out-of-place transformation
-
-

-

-

- Of course, you do not need to call these two method. You can always use - the simpler do_step( - sys , - inout , - t , - dt ), - but sometimes the derivative of the state is needed externally to do some - external computations or to perform some statistical analysis. -

-

- A special class of the explicit steppers are the FSAL (first-same-as-last) - steppers, where the last evaluation of the system function is also the - first evaluation of the following step. For such steppers the do_step method are slightly different: -

-

- do_step( - sys , - inout , - dxdtinout , - t , - dt ) -

-

- do_step( - sys , - in , - dxdtin , - out , - dxdtout , - t , - dt ) -

-

- This method takes the derivative at time t - and also also stores the derivative at time t+dt. - Calling these functions subsequently iterating along the solution one saves - one function call by passing the result for dxdt into the next function - call. However, when using FSAL steppers without supplying derivatives: -

-

- do_step( - sys , - inout , - t , - dt ) -

-

- the stepper internally satisfies the FSAL property which means it remembers - the last dxdt and uses - it for the next step. An example for a FSAL stepper is the Runge-Kutta-Dopri5 - stepper. The FSAL-trick is sometimes also referred as the Fehlberg-Trick. - An example how the FSAL steppers can be used is -

-

-

-
runge_kutta_dopri5< state_type > rk;
-rk.do_step( sys1 , in , t , out , dt );
-rk.do_step( sys2 , in , t , out , dt );         // DONT do this, sys1 is assumed
-
-rk.do_step( sys2 , in2 , t , out , dt );
-rk.do_step( sys2 , in3 , t , out , dt );        // DONT do this, in2 is assumed
-
-rk.do_step( sys1 , inout , dxdtinout , t , dt );
-rk.do_step( sys2 , inout , dxdtinout , t , dt );           // Ok, internal derivative is not used, dxdtinout is updated
-
-rk.do_step( sys1 , in , dxdtin , t , out , dxdtout , dt );
-rk.do_step( sys2 , in , dxdtin , t , out , dxdtout , dt ); // Ok, internal derivative is not used
-
-

-

-
- - - - - -
[Caution]Caution

- The FSAL-steppers save the derivative at time t+dt - internally if they are called via do_step( sys , in , out , t , dt ). The first call of do_step - will initialize dxdt - and for all following calls it is assumed that the same system and the - same state are used. If you use the FSAL stepper within the integrate - functions this is taken care of automatically. See the Using - steppers section for more details or look into the table below - to see which stepper have an internal state. -

-
-
-

-Symplectic - solvers -

-

- As mentioned above symplectic solvers are used for Hamiltonian systems. - Symplectic solvers conserve the phase space volume exactly and if the Hamiltonian - system is energy conservative they also conserve the energy approximately. - A special class of symplectic systems are separable systems which can be - written in the form dqdt/dt = f1(p), dpdt/dt - = f2(q), where (q,p) are the state of system. - The space of (q,p) is sometimes referred as the phase - space and q and p are said the - be the phase space variables. Symplectic systems in this special form occur - widely in nature. For example the complete classical mechanics as written - down by Newton, Lagrange and Hamilton can be formulated in this framework. - Of course, the separability of the system depends on the specific choice - of coordinates. -

-

- Symplectic systems can be solved by odeint by means of the symplectic_euler - stepper and a symplectic Runge-Kutta-Nystrom method of sixth-order. These - steppers assume that the system is autonomous, hence the time will not - explicitly occur. Further they fulfill in principle the default Stepper - concept, but they expect the system to be a pair of callable objects. The - first entry of this pair calculates f1(p) while the - second calculates f2(q). The syntax is sys.first(p,dqdt) and sys.second(q,dpdt), - where the first and second part can be again simple C-functions of functors. - An example is the harmonic oscillator: -

-

-

-
typedef boost::array< double , 1 > vector_type;
-
-
-struct harm_osc_f1
-{
-    void operator()( const vector_type &p , vector_type &dqdt )
-    {
-        dqdt[0] = p[0];
-    }
-};
-
-struct harm_osc_f2
-{
-    void operator()( const vector_type &q , vector_type &dpdt )
-    {
-        dpdt[0] = -q[0];
-    }
-};
-
-

-

-

- The state of such an ODE consist now also of two parts, the part for q - (also called the coordinates) and the part for p (the momenta). The full - example for the harmonic oscillator is now: -

-

-

-
pair< vector_type , vector_type > x;
-x.first[0] = 1.0; x.second[0] = 0.0;
-symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
-rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , x , t , dt );
-
-

-

-

- If you like to represent the system with one class you can easily bind - two public method: -

-

-

-
struct harm_osc
-{
-    void f1( const vector_type &p , vector_type &dqdt ) const
-    {
-        dqdt[0] = p[0];
-    }
-
-    void f2( const vector_type &q , vector_type &dpdt ) const
-    {
-        dpdt[0] = -q[0];
-    }
-};
-
-

-

-

-

-
harm_osc h;
-rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) ,
-        x , t , dt );
-
-

-

-

- Many Hamiltonian system can be written as dq/dt=p, - dp/dt=f(q) which is computationally much easier than - the full separable system. Very often, it is also possible to transform - the original equations of motion to bring the system in this simplified - form. This kind of system can be used in the symplectic solvers, by simply - passing f(p) to the do_step - method, again f(p) will be represented by a simple - C-function or a functor. Here, the above example of the harmonic oscillator - can be written as -

-

-

-
pair< vector_type , vector_type > x;
-x.first[0] = 1.0; x.second[0] = 0.0;
-symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
-rkn.do_step( harm_osc_f1() , x , t , dt );
-
-

-

-

- In this example the function harm_osc_f1 - is exactly the same function as in the above examples. -

-

- Note, that the state of the ODE must not be constructed explicitly via - pair< - vector_type , - vector_type > - x. One can also use a combination - of make_pair and ref. Furthermore, a convenience version - of do_step exists which - takes q and p without combining them into a pair: -

-

-

-
rkn.do_step( harm_osc_f1() , make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt );
-rkn.do_step( harm_osc_f1() , q , p , t , dt );
-rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , q , p , t , dt );
-
-

-

-
-
-

-Implicit - solvers -

-
- - - - - -
[Caution]Caution

- This section is not up-to-date. -

-

- For some kind of systems the stability properties of the classical Runge-Kutta - are not sufficient, especially if the system is said to be stiff. A stiff - system possesses two or more time scales of very different order. Solvers - for stiff systems are usually implicit, meaning that they solve equations - like x(t+dt) = x(t) + dt * f(x(t+1)). This particular - scheme is the implicit Euler method. Implicit methods usually solve the - system of equations by a root finding algorithm like the Newton method - and therefore need to know the Jacobian of the system J​ij = df​i / - dx​j. -

-

- For implicit solvers the system is again a pair, where the first component - computes f(x,t) and the second the Jacobian. The syntax - is sys.first( x , dxdt , t ) and - sys.second( x , J , t ). - For the implicit solver the state_type - is ublas::vector and the Jacobian is represented - by ublas::matrix. -

-
- - - - - -
[Important]Important

- Implicit solvers do only work with ublas::vector as state type. At the - moment, no other state types are supported. -

-
-
-

-Multistep - methods -

-

- Another large class of solvers are multi-step method. They save a small - part of the history of the solution and compute the next step with the - help of this history. Since multi-step methods know a part of their history - they do not need to compute the system function very often, usually it - is only computed once. This makes multi-step methods preferable if a call - of the system function is expensive. Examples are ODEs defined on networks, - where the computation of the interaction is usually where expensive (and - might be of order O(N^2)). -

-

- Multi-step methods differ from the normal steppers. They safe a part of - their history and this part has to be explicitly calculated and initialized. - In the following example an Adams-Bashforth-stepper with a history of 5 - steps is instantiated and initialized; -

-

-

-
adams_bashforth_moulton< 5 , state_type > abm;
-abm.initialize( sys , inout , t , dt );
-abm.do_step( sys , inout , t , dt );
-
-

-

-

- The initialization uses a fourth-order Runge-Kutta stepper and after the - call of initialize the - state of inout has changed - to the current state, such that it can be immediately used by passing it - to following calls of do_step. - Of course, you can also use you own steppers to initialize the internal - state of the Adams-Bashforth-Stepper: -

-

-

-
abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt );
-
-

-

-

- Many multi-step methods are also explicit steppers, hence the parameter - of do_step method do not - differ from the explicit steppers. -

-
- - - - - -
[Caution]Caution

- The multi-step methods have some internal variables which depend on the - explicit solution. Hence after any external changes of your state (e.g. - size) or system the initialize function has to be called again to adjust - the internal state of the stepper. Of course, if you use the integrate - functions this will be taken into account. See the Using - steppers section for more details. -

-
-
-

-Controlled - steppers -

-

- Many of the above introduced steppers possess the possibility to use adaptive - step-size control. Adaptive step size integration works in principle as - follows: -

-
    -
  1. - The error of one step is calculated. This is usually done by performing - two steps with different orders. The difference between these two steps - is then used as a measure for the error. Stepper which can calculate - the error are Error - Stepper and they form an own class with an separate concept. -
    2. -
  2. - This error is compared against some predefined error tolerances. Are - the tolerance violated the step is reject and the step-size is decreases. - Otherwise the step is accepted and possibly the step-size is increased. -
    3. -
-

- The class of controlled steppers has their own concept in odeint - the - Controlled - Stepper concept. They are usually constructed from the underlying - error steppers. An example is the controller for the explicit Runge-Kutta - steppers. The Runge-Kutta steppers enter the controller as a template argument. - Additionally one can pass the Runge-Kutta stepper to the constructor, but - this step is not necessary; the stepper is default-constructed if possible. -

-

- Different step size controlling mechanism exist. They all have in common - that they somehow compare predefined error tolerance against the error - and that they might reject or accept a step. If a step is rejected the - step size is usually decreased and the step is made again with the reduced - step size. This procedure is repeated until the step is accepted. This - algorithm is implemented in the integration functions. -

-

- A classical way to decide whether a step is rejected or accepted is to - calculate -

-

- val = || | err​i | / ( ε​abs + ε​rel * ( a​x | x​i | + a​dxdt | | dxdt​i | )|| - -

-

- ε​abs and ε​rel are the absolute - and the relative error tolerances, and || x || is - a norm, typically ||x||=(Σ​i x​i2)1/2 or the maximum norm. - The step is rejected if val is greater then 1, otherwise - it is accepted. For details of the used norms and error tolerance see the - table below. -

-

- For the controlled_runge_kutta - stepper the new step size is then calculated via -

-

- val > 1 : dt​new = dt​current max( 0.9 pow( val , -1 / ( O​E - 1 - ) ) , 0.2 ) -

-

- val < 0.5 : dt​new = dt​current min( 0.9 pow( val , -1 / O​S ) , - 5 ) -

-

- else : dt​new = dt​current -

-

- Here, O​S and O​E are the order - of the stepper and the error stepper. These formulas also contain some - safety factors, avoiding that the step size is reduced or increased to - much. For details of the implementations of the controlled steppers in - odeint see the table below. -

-
-

Table 1.6. Adaptive step size algorithms

-
------ - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Stepper -

- 		-

- Tolerance formula -

- 		-

- Norm -

- 		-

- Step size adaption -

-
-

- controlled_runge_kutta -

- 		-

- val = || | err​i | / ( ε​abs + ε​rel * ( a​x | x​i | + a​dxdt | | - dxdt​i | )|| -

- 		-

- ||x|| = max( x​i ) -

- 		-

- val > 1 : dt​new = dt​current max( 0.9 pow( val , -1 - / ( O​E - 1 ) ) , 0.2 ) -

-

- val < 0.5 : dt​new = dt​current min( 0.9 pow( val , - -1 / O​S ) , 5 ) -

-

- else : dt​new = dt​current -

-
-

- rosenbrock4_controller -

- 		-

- val = || err​i / ( ε​abs + ε​rel max( | x​i | , | xold​i | ) ) - || -

- 		-

- ||x||=(Σ​i x​i2)1/2 -

- 		-

- fac = max( 1 / 6 , min( 5 , pow( val , 1 / 4 ) / 0.9 - ) -

-

- fac2 = max( 1 / 6 , min( 5 , dt​old / dt​current pow( val2 / - val​old , 1 / 4 ) / 0.9 ) -

-

- val > 1 : dt​new = dt​current / fac -

-

- val < 1 : dt​new = dt​current / max( fac , fac2 ) -

-
-

- bulirsch_stoer -

- 		-

- tol=1/2 -

- 		-

- - -

- 		-

- dt​new = dt​old1/a -

-
-
-

- safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0 -

-

- value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , - 0.25 ) / safe; fac_pred = std::max( fac2 , std::min( fac1 , fac_pred ) - ); fac = std::max( fac , fac_pred ); dt_new = dt / fac; -

-

- fac = max( fac2 , min( fac1 , pow( err , 0.25 ) / safe ) ) -

-

- To ease to generation of the controlled stepper, generation functions exist - which take the absolute and relative error tolerances and a predefined - error stepper and construct from this knowledge an appropriate controlled - stepper. The generation functions are explained in detail in Generation - functions. -

-
-
-

-Dense - output steppers -

-

- A fourth class of stepper exists which are the so called dense output steppers. - Dense-output steppers might take larger steps and interpolate the solution - between two consecutive points. This interpolated points have usually the - same order as the order of the stepper. Dense-output steppers are often - composite stepper which take the underlying method as a template parameter. - An example is the dense_output_runge_kutta - stepper which takes a Runge-Kutta stepper with dense-output facilities - as argument. Not all Runge-Kutta steppers provide dense-output calculation; - at the moment only the Dormand-Prince 5 stepper provides dense output. - An example is -

-

-

-
dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense;
-dense.initialize( in , t , dt );
-pair< double , double > times = dense.do_step( sys );
-
-

-

-

- Dense output stepper have their own concept. The main difference to usual - steppers is that they manage the state and time internally. If you call - do_step, only the ODE is - passed as argument. Furthermore do_step - return the last time interval: t - and t+dt, hence you can interpolate the solution - between these two times points. Another difference is that they must be - initialized with initialize, - otherwise the internal state of the stepper is default constructed which - might produce funny errors or bugs. -

-

- The construction of the dense output stepper looks a little bit nasty, - since in the case of the dense_output_runge_kutta - stepper a controlled stepper and an error stepper have to be nested. To - simplify the generation of the dense output stepper generation functions - exist: -

-

-

-
boost::numeric::odeint::result_of::make_dense_output< runge_kutta_dopri5< state_type > >::type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
-
-

-

-

- Of course, this statement is also lengthy; it demonstrates how make_dense_output can be used with the - result_of protocol. The - parameters to make_dense_output - are the absolute error tolerance, the relative error tolerance and the - stepper. This explicitly assumes that the underlying stepper is a controlled - stepper and that this stepper has an absolute and a relative error tolerance. - For details about the generation functions see Generation - functions. The generation functions have been designed for easy - use with the integrate functions: -

-

-

-
integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt );
-
-

-

-
-
-

-Using - steppers -

-

- This section contains some general information about the usage of the steppers - in odeint. -

-

- Steppers are copied by value -

-

- The stepper in odeint are always copied by values. They are copied for - the creation of the controlled steppers or the dense output steppers as - well as in the integrate functions. -

-

- Steppers might have a internal state -

-
- - - - - -
[Caution]Caution

- Some of the features described in this section are not yet implemented -

-

- Some steppers require to store some information about the state of the - ODE between two steps. Examples are the multi-step methods which store - a part of the solution during the evolution of the ODE, or the FSAL steppers - which store the last derivative at time t+dt, to be - used in the next step. In both cases the steppers expect that consecutive - calls of do_step are from - the same solution and the same ODE. In this case it is absolutely necessary - that you call do_step with - the same system function and the same state, see also the examples for - the FSAL steppers above. -

-

- Stepper with an internal state support two additional methods: reset which resets the state and initialize which initializes the internal - state. The parameters of initialize - depend on the specific stepper. For example the Adams-Bashforth-Moulton - stepper provides two initialize methods: initialize( system , inout , t , dt ) which initializes the internal states - with the help of the Runge-Kutta 4 stepper, and initialize( stepper , system , inout , t , dt ) which initializes with the help of stepper. For the case of the FSAL steppers, - initialize is initialize( - sys , - in , - t ) - which simply calculates the r.h.s. of the ODE and assigns its value to - the internal derivative. -

-

- All these steppers have in common, that they initially fill their internal - state by themselves. Hence you are not required to call initialize. See - how this works for the Adams-Bashforth-Moulton stepper: in the example - we instantiate a fourth order Adams-Bashforth-Moulton stepper, meaning - that it will store 4 internal derivatives of the solution at times (t-dt,t-2*dt,t-3*dt,t-4*dt). -

-

-

-
adams_bashforth_moulton< 4 , state_type > stepper;
-stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the first internal state
-                                       // the internal array is now [x(t-dt)]
-
-stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the second internal state
-                                       // the internal state array is now [x(t-dt), x(t-2*dt)]
-
-stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the third internal state
-                                       // the internal state array is now [x(t-dt), x(t-2*dt), x(t-3*dt)]
-
-stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the fourth internal state
-                                       // the internal state array is now [x(t-dt), x(t-2*dt), x(t-3*dt), x(t-4*dt)]
-
-stepper.do_step( sys , x , t , dt );   // make one step with Adam-Bashforth-Moulton, the internal array of states is now rotated
-
-

-

-

- In the stepper table at the bottom of this page one can see which stepper - have an internal state and hence provide the reset - and initialize methods. -

-

- Stepper might be resizable -

-

- Nearly all steppers in odeint need to store some intermediate results of - the type state_type or - deriv_type. To do so odeint - need some memory management for the internal temporaries. As this memory - management is typically related to adjusting the size of vector-like types, - it is called resizing in odeint. So, most steppers in odeint provide an - additional template parameter which controls the size adjustment of the - internal variables - the resizer. In detail odeint provides three policy - classes (resizers) always_resizer, - initially_resizer, and - never_resizer. Furthermore, - all stepper have a method adjust_size - which takes a parameter representing a state type and which manually adjusts - the size of the internal variables matching the size of the given instance. - Before performing the actual resizing odeint always checks if the sizes - of the state and the internal variable differ and only resizes if they - are different. -

-
- - - - - -
[Note]Note

- You only have to worry about memory allocation when using dynamically - sized vector types. If your state type is heap allocated, like boost::array, no memory allocation is required - whatsoever. -

-

- By default the resizing parameter is initially_resizer, - meaning that the first call to do_step - performs the resizing, hence memory allocation. Of course, if you have - changed the size of your system and your state you have to call adjust_size by hand in this case. The - second resizer is the always_resizer - which tries to resize the internal variables at every call of do_step. Typical use cases for this kind - of resizer are self expanding lattices like shown in the tutorial or partial - differential equations with an adaptive grid. Here, no calls of adjust_size are required, the steppers - manage everything themselve. The third class of resizer is the never_resizer which means that the internal - variables are never adjusted automatically and always have to be adjusted - by hand . -

-

- There is a second mechanism which influences the resizing and which controls - if a state type is at least resizeable - a meta-function is_resizeable. This meta-function returns - a static boolean value if any type is resizable. For example it will return - true for std::vector< T > but false - for boost::array< T >. - By default and for unknown types is_resizeable - returns false, so if you have - your own type you need to specialize this meta-function. For more details - on the resizing mechanism see the section Adapt - you own state types. -

-

- Which steppers should be used in which situation -

-

- odeint provides a quite large number of different steppers such that the - user is left with the question of which stepper fits his needs. Our personal - recommendations are: -

-
    -
  • - runge_kutta_dopri5 - is maybe the best default stepper. It has step size control as well - as dense-output functionality. Simple create a dense-output stepper - by make_dense_output( 1.0e-6 , 1.0e-5 , runge_kutta_dopri5< state_type - >() ). -
    • -
  • - runge_kutta4 is a good - stepper for constant step sizes. It is widely used and very well known. - If you need to create artificial time series this stepper should be - the first choice. -
    • -
  • - 'runge_kutta_fehlberg78' is similar to the 'runge_kutta4' with the - advantage that it has higher precision. It can also be used with step - size control. -
    • -
  • - adams_bashforth_moulton - is very well suited for ODEs where the r.h.s. is expensive (in terms - of computation time). It will calculate the system function only once - during each step. -
    • -
-
-
-

-Stepper - overview -

-
-

Table 1.7. Stepper Algorithms

-
----------- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Algorithm -

- 		-

- Class -

- 		-

- Concept -

- 		-

- System Concept -

- 		-

- Order -

- 		-

- Error Estimation -

- 		-

- Dense Output -

- 		-

- Internal state -

- 		-

- Remarks -

-
-

- Explicit Euler -

- 		-

- euler -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- 1 -

- 		-

- No -

- 		-

- Yes -

- 		-

- No -

- 		-

- Very simple, only for demonstrating purpose -

-
-

- Modified Midpoint -

- 		-

- modified_midpoint -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable (2) -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Used in Bulirsch-Stoer implementation -

-
-

- Runge-Kutta 4 -

- 		-

- runge_kutta4 -

- 		-

- Stepper -

- 		-

- System -

- 		-

- 4 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- The classical Runge Kutta scheme, good general scheme without - error control -

-
-

- Cash-Karp -

- 		-

- runge_kutta_cash_karp54 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 5 -

- 		-

- Yes (4) -

- 		-

- No -

- 		-

- No -

- 		-

- Good general scheme with error estimation, to be used in controlled_error_stepper -

-
-

- Dormand-Prince 5 -

- 		-

- runge_kutta_dopri5 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 5 -

- 		-

- Yes (4) -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- Standard method with error control and dense output, to be used - in controlled_error_stepper and in dense_output_controlled_explicit_fsal. -

-
-

- Fehlberg 78 -

- 		-

- runge_kutta_fehlberg78 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 8 -

- 		-

- Yes (7) -

- 		-

- No -

- 		-

- No -

- 		-

- Good high order method with error estimation, to be used in controlled_error_stepper. -

-
-

- Adams Bashforth -

- 		-

- adams_bashforth -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Multistep method -

-
-

- Adams Moulton -

- 		-

- adams_moulton -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Multistep method -

-
-

- Adams Bashforth Moulton -

- 		-

- adams_bashforth_moulton -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Combined multistep method -

-
-

- Controlled Runge Kutta -

- 		-

- controlled_runge_kutta -

- 		-

- Controlled - Stepper -

- 		-

- System -

- 		-

- depends -

- 		-

- Yes -

- 		-

- No -

- 		-

- depends -

- 		-

- Error control for Error - Stepper. Requires an Error - Stepper from above. Order depends on the given ErrorStepper -

-
-

- Dense Output Runge Kutta -

- 		-

- dense_output_runge_kutta -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- depends -

- 		-

- No -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- Dense output for Stepper - and Error - Stepper from above if they provide dense output functionality - (like euler and - runge_kutta_dopri5). - Order depends on the given stepper. -

-
-

- Bulirsch-Stoer -

- 		-

- bulirsch_stoer -

- 		-

- Controlled - Stepper -

- 		-

- System -

- 		-

- variable -

- 		-

- Yes -

- 		-

- No -

- 		-

- No -

- 		-

- Stepper with step size and order control. Very good if high precision - is required. -

-
-

- Bulirsch-Stoer Dense Output -

- 		-

- bulirsch_stoer_dense_out -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- variable -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Stepper with step size and order control as well as dense output. - Very good if high precision and dense output is required. -

-
-

- Implicit Euler -

- 		-

- implicit_euler -

- 		-

- Stepper -

- 		-

- Implicit - System -

- 		-

- 1 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Basic implicit routine. Requires the Jacobian. Works only with - Boost.UBlas - vectors as state types. -

-
-

- Rosenbrock 4 -

- 		-

- rosenbrock4 -

- 		-

- Error - Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Good for stiff systems. Works only with Boost.UBlas - vectors as state types. -

-
-

- Controlled Rosenbrock 4 -

- 		-

- rosenbrock4_controller -

- 		-

- Controlled - Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Rosenbrock 4 with error control. Works only with Boost.UBlas - vectors as state types. -

-
-

- Dense Output Rosenbrock 4 -

- 		-

- rosenbrock4_dense_output -

- 		-

- Dense - Output Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Controlled Rosenbrock 4 with dense output. Works only with Boost.UBlas - vectors as state types. -

-
-

- Symplectic Euler -

- 		-

- symplectic_euler -

- 		-

- Stepper -

- 		-

- Symplectic - System Simple - Symplectic System -

- 		-

- 1 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Basic symplectic solver for separable Hamiltonian system -

-
-

- Symplectic RKN McLachlan -

- 		-

- symplectic_rkn_sb3a_mclachlan -

- 		-

- Stepper -

- 		-

- Symplectic - System Simple - Symplectic System -

- 		-

- 6 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Symplectic solver for separable Hamiltonian system with order - 6 -

-
-
-
-
-
-

-Write - own steppers -

-

- Finally, one can also write new steppers which are fully compatible with - odeint. They only have to fulfill one or several of the stepper Concepts - of odeint. -

-

- We will illustrate how to write your own stepper with the example of the - stochastic Euler method. This method is suited to solve stochastic differential - equations (SDEs). A SDE has the form -

-

- dx/dt = f(x) + g(x) ξ(t) -

-

- where ξ is Gaussian white noise with zero mean and - a standard deviation σ(t). f(x) - is said to be the deterministic part while g(x) ξ is - the noisy part. In case g(x) is independent of x - the SDE is said to have additive noise. It is not possible to solve SDE - with the classical solvers for ODEs since the noisy part of the SDE has - to be scaled differently then the deterministic part with respect to the - time step. But there exist many solvers for SDEs. A classical and easy - method is the stochastic Euler solver. It works by iterating -

-

- x(t+Δ t) = x(t) + Δ t f(x(t)) + Δ t1/2 g(x) ξ(t) -

-

- where ξ(t) is an independent normal distributed random variable. -

-

- Now we will implement this method. We will call the stepper stochastic_euler. It models the Stepper concept. - For simplicity, we fix the state type to be an array< double , N > The class definition looks like -

-

-

-
template< size_t N > class stochastic_euler
-{
-public:
-
-    typedef boost::array< double , N > state_type;
-    typedef boost::array< double , N > deriv_type;
-    typedef double value_type;
-    typedef double time_type;
-    typedef unsigned short order_type;
-    typedef boost::numeric::odeint::stepper_tag stepper_category;
-
-    static order_type order( void ) { return 1; }
-
-    // ...
-};
-
-

-

-

- The types are needed in order to fulfill the stepper concept. As internal - state and deriv type we use simple arrays in the stochastic Euler, they - are needed for the temporaries. The stepper has the order one which is - returned from the order() function. -

-

- The system functions needs to calculate the deterministic and the stochastic - part of our stochastic differential equation. So it might be suitable that - the system function itself is a pair of functions, one for computing the - deterministic and one for the stochastic part. The first element of the - pair simply computes the deterministic part while the second the stochastic - one. Then, the second part also needs to calculate the random numbers in - order to simulate the stochastic process. We can now implement the do_step method -

-

-

-
template< size_t N > class stochastic_euler
-{
-public:
-
-    // ...
-
-    template< class System >
-    void do_step( System system , state_type &x , time_type t , time_type dt ) const
-    {
-        deriv_type det , stoch ;
-        system.first( x , det );
-        system.second( x , stoch );
-        for( size_t i=0 ; i<x.size() ; ++i )
-            x[i] += dt * det[i] + sqrt( dt ) * stoch[i];
-    }
-};
-
-

-

-

- This is all. It is quite simple and the stochastic Euler stepper implement - here is quite general. Of course it can be enhanced, for example -

-
    -
  • - use of operations and algebras as well as the resizing mechanism for - maximal flexibility and portability -
    • -
  • - use of boost::ref for the system functions -
    • -
  • - use of boost::range for the state type in the - do_step method -
    • -
  • - ... -
    • -
-

- Now, lets look how we use the new stepper. A nice example is the Ornstein-Uhlenbeck - process. It consists of a simple Brownian motion overlapped with an relaxation - process. Its SDE reads -

-

- dx/dt = - x + ξ -

-

- where ξ is Gaussian white noise with standard deviation σ. - Implementing the Ornstein-Uhlenbeck process is quite simple. We need two - functions or functors - one for the deterministic and one for the stochastic - part: -

-

-

-
const static size_t N = 1;
-typedef boost::array< double , N > state_type;
-
-struct ornstein_det
-{
-    void operator()( const state_type &x , state_type &dxdt ) const
-    {
-        dxdt[0] = -x[0];
-    }
-};
-
-struct ornstein_stoch
-{
-    boost::mt19937 m_rng;
-    boost::normal_distribution<> m_dist;
-
-    ornstein_stoch( double sigma ) : m_rng() , m_dist( 0.0 , sigma ) { }
-
-    void operator()( const state_type &x , state_type &dxdt )
-    {
-        dxdt[0] = m_dist( m_rng );
-    }
-};
-
-

-

-

- In the stochastic part we have used the Mersenne twister for the random - number generation and a Gaussian white noise generator normal_distribution - with standard deviation σ. Now, we can use the stochastic - Euler stepper with the integrate functions: -

-

-

-
double dt = 0.1;
-state_type x = {{ 1.0 }};
-integrate_const( stochastic_euler< N >() , make_pair( ornstein_det() , ornstein_stoch( 1.0 ) ) ,
-        x , 0.0 , 10.0 , dt , streaming_observer() );
-
-

-

-

- Note, how we have used the make_pair - function for the generation of the system function. -

-
-
-

-Writing - own Runge-Kutta steppers -

-

- odeint provides a C++ template meta-algorithm for constructing arbitrary - Runge-Kutta schemes [1]. Some schemes are predefined in odeint, for example the classical - Runge-Kutta of fourth order, or the Runge-Kutta-Cash-Karp 54 and the Runge-Kutta-Fehlberg - 78 method. Of course, you can use this meta algorithm to construct you - own solvers. This has the advantage that you can make full use of odeint's - algebra and operation system. -

-

- Consider for example the method of Heun, defined by the following Butcher - tableau: -

-
c1 = 0
-
-c2 = 1/3, a21 = 1/3
-
-c3 = 2/3, a31 =  0 , a32 = 2/3
-
-          b1  = 1/4, b2  = 0  , b3 = 3/4
-
-

- Implementing this method is very easy. First you have to define the constants: -

-

-

-
template< class Value = double >
-struct heun_a1 : boost::array< Value , 1 > {
-    heun_a1( void )
-    {
-        (*this)[0] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
-    }
-};
-
-template< class Value = double >
-struct heun_a2 : boost::array< Value , 2 >
-{
-    heun_a2( void )
-    {
-        (*this)[0] = static_cast< Value >( 0 );
-        (*this)[1] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
-    }
-};
-
-
-template< class Value = double >
-struct heun_b : boost::array< Value , 3 >
-{
-    heun_b( void )
-    {
-        (*this)[0] = static_cast<Value>( 1 ) / static_cast<Value>( 4 );
-        (*this)[1] = static_cast<Value>( 0 );
-        (*this)[2] = static_cast<Value>( 3 ) / static_cast<Value>( 4 );
-    }
-};
-
-template< class Value = double >
-struct heun_c : boost::array< Value , 3 >
-{
-    heun_c( void )
-    {
-        (*this)[0] = static_cast< Value >( 0 );
-        (*this)[1] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
-        (*this)[2] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
-    }
-};
-
-

-

-

- Ok, you are right, this looks rather ugly. Nevertheless, packing all parameters - into a templatized class which is not immediately evaluated has the advantage - that you can change the value_type - of your stepper to any type you like - presumably arbitrary precision types. - Of course, one could instantiate the coefficients directly -

-

-

-
const boost::array< double , 1 > heun_a1 = {{ 1.0 / 3.0 }};
-const boost::array< double , 2 > heun_a2 = {{ 0.0 , 2.0 / 3.0 }};
-const boost::array< double , 3 > heun_b = {{ 1.0 / 4.0 , 0.0 , 3.0 / 4.0 }};
-const boost::array< double , 3 > heun_c = {{ 0.0 , 1.0 / 3.0 , 2.0 / 3.0 }};
-
-

-

-

- But then you are nailed down to use doubles. -

-

- Next, you need to define your stepper, note that the Heun method has 3 - stages and produces approximations of order 3: -

-

-

-
template<
-    class State ,
-    class Value = double ,
-    class Deriv = State ,
-    class Time = Value ,
-    class Algebra = boost::numeric::odeint::range_algebra ,
-    class Operations = boost::numeric::odeint::default_operations ,
-    class Resizer = boost::numeric::odeint::initially_resizer
->
-class heun : public
-boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
-                                             Algebra , Operations , Resizer >
-{
-
-public:
-
-    typedef boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
-                                                         Algebra , Operations , Resizer > stepper_base_type;
-
-    typedef typename stepper_base_type::state_type state_type;
-    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
-    typedef typename stepper_base_type::value_type value_type;
-    typedef typename stepper_base_type::deriv_type deriv_type;
-    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
-    typedef typename stepper_base_type::time_type time_type;
-    typedef typename stepper_base_type::algebra_type algebra_type;
-    typedef typename stepper_base_type::operations_type operations_type;
-    typedef typename stepper_base_type::resizer_type resizer_type;
-    typedef typename stepper_base_type::stepper_type stepper_type;
-
-    heun( const algebra_type &algebra = algebra_type() )
-    : stepper_base_type(
-            fusion::make_vector(
-                heun_a1<Value>() ,
-                heun_a2<Value>() ) ,
-            heun_b<Value>() , heun_c<Value>() , algebra )
-    { }
-};
-
-

-

-

- That's it. Now, we have a new stepper method and we can use it, for example - with the Lorenz system: -

-

-

-
typedef boost::array< double , 3 > state_type;
-heun< state_type > h;
-state_type x = {{ 10.0 , 10.0 , 10.0 }};
-
-integrate_const( h , lorenz() , x , 0.0 , 100.0 , 0.01 ,
-                 streaming_observer( std::cout ) );
-
-

-

-
-
-
-

[1] - M. Mulansky, K. Ahnert, Template-Metaprogramming applied to numerical - problems, arxiv:1110.3233 -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/odeint_in_detail/using_boost__range.html b/doc/boost_numeric_odeint/odeint_in_detail/using_boost__range.html deleted file mode 100644 index 99e19f34..00000000 --- a/doc/boost_numeric_odeint/odeint_in_detail/using_boost__range.html +++ /dev/null @@ -1,246 +0,0 @@ - - - -Using boost::range - - - - - - - - -
-PrevUpHomeNext -
-
-

-Using - boost::range -

-

- Most steppers in odeint also accept the state give as a range. A range is - sequence of values modeled by a range concept. See Boost.Range - for an overview over existing concepts and examples of ranges. This means - that the state_type of the - stepper must not necessarily be used to call the do_step - method. -

-

- One use-case for Boost.Range - in odeint has been shown in Chaotic - System where the state consists of two parts: one for the original - system and one for the perturbations. The ranges are used to initialize (solve) - only the system part where the perturbation part is not touched, that is - a range consisting only of the system part is used. After that the complete - state including the perturbations is solved. -

-

- Another use case is a system consisting of coupled units where you want to - initialize each unit separately with the ODE of the uncoupled unit. An example - is a chain of coupled van-der-Pol-oscillators which are initialized uniformly - from the uncoupled van-der-Pol-oscillator. Then you can use Boost.Range - to solve only one individual oscillator in the chain. -

-

- In short, you can Boost.Range - to use one state within two system functions which expect states with different - sizes. -

-

- An example was given in the Chaotic - System tutorial. Using Boost.Range usually means that your system - function needs to adapt to the iterators of Boost.Range. That is, your function - is called with a range and you need to get the iterators from that range. - This can easily be done. You have to implement your system as a class or - a struct and you have to templatize the operator(). Then you can use the range_iterator-meta - function and boost::begin and boost::end to - obtain the iterators of your range: -

-

-

-
class sys
-{
-    template< class State , class Deriv >
-    void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
-    {
-         typename boost::range_iterator< const State >::type x = boost::begin( x_ );
-         typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
-
-         // fill dxdt
-    }
-};
-
-

-

-

- If your range is a random access-range you can also apply the bracket operator - to the iterator to access the elements in the range: -

-
class sys
-{
-    template< class State , class Deriv >
-    void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
-    {
-         typename boost::range_iterator< const State >::type x = boost::begin( x_ );
-         typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
-
-         dxdt[0] = f1( x[0] , x[1] );
-         dxdt[1] = f2( x[0] , x[1] );
-    }
-};
-
-

-

-

- The following two tables show which steppers and which algebras are compatible - with Boost.Range. -

-
-

Table 1.10. Steppers supporting Boost.Range

-
- - - - - - - - - - - - - - - - - - - - - - - -
-

- Stepper -

-
-

- adams_bashforth_moulton -

-
-

- bulirsch_stoer_dense_out -

-
-

- bulirsch_stoer -

-
-

- controlled_runge_kutta -

-
-

- dense_output_runge_kutta -

-
-

- euler -

-
-

- explicit_error_generic_rk -

-
-

- explicit_generic_rk -

-
-

- rosenbrock4_controller -

-
-

- rosenbrock4_dense_output -

-
-

- rosenbrock4 -

-
-

- runge_kutta4_classic -

-
-

- runge_kutta4 -

-
-

- runge_kutta_cash_karp54_classic -

-
-

- runge_kutta_cash_karp54 -

-
-

- runge_kutta_dopri5 -

-
-

- runge_kutta_fehlberg78 -

-
-

- symplectic_euler -

-
-

- symplectic_rkn_sb3a_mclachlan -

-
-
-
-

Table 1.11. Algebras supporting Boost.Range

-
- - - - - - -
-

- algebra -

-
-

- range_algebra -

-
-

- thrust_algebra -

-
-
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/odeint_in_detail/using_boost__ref.html b/doc/boost_numeric_odeint/odeint_in_detail/using_boost__ref.html deleted file mode 100644 index 157c1aab..00000000 --- a/doc/boost_numeric_odeint/odeint_in_detail/using_boost__ref.html +++ /dev/null @@ -1,79 +0,0 @@ - - - -Using boost::ref - - - - - - - - -
-PrevUpHomeNext -
-
-

-Using - boost::ref -

-

- In odeint all system functions and observers are passed by value. For example, - if you call a do_step method - of a particular stepper or the integration functions, your system and your - stepper will be passed by value: -

-

-

-
rk4.do_step( sys , x , t , dt );   // pass sys by value
-
-

-

-

- This behavior is suitable for most systems, especially if your system does - not contain any data or only a few parameters. However, in some cases you - might contain some large amount of data with you system function and passing - them by value is not desired since the data would be copied. -

-

- In such cases you can easily use boost::ref (and - its relative boost::cref) which passes its argument by reference - (or constant reference). odeint will unpack the arguments and no copying - at all of your system object will take place: -

-

-

-
rk4.do_step( boost::ref( sys ) , x , t , dt );   // pass sys as references
-
-

-

-

- The same mechanism can be used for the observers in the integrate functions. -

-
- - - - - -
[Tip]Tip

- If you are using C++11 you can also use std::ref - and std::cref -

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial.html b/doc/boost_numeric_odeint/tutorial.html deleted file mode 100644 index 4bbd6a0b..00000000 --- a/doc/boost_numeric_odeint/tutorial.html +++ /dev/null @@ -1,51 +0,0 @@ - - - -Tutorial - - - - - - - - -
-PrevUpHomeNext -
-
-

-Tutorial -

-
-
Harmonic - oscillator
-
Solar system
-
Chaotic - systems and Lyapunov exponents
-
Stiff systems
-
Special - topics
-
Using - Cuda and Thrust
-
Using - OpenCL and VexCL
-
All examples
-
References
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/all_examples.html b/doc/boost_numeric_odeint/tutorial/all_examples.html deleted file mode 100644 index f013bee5..00000000 --- a/doc/boost_numeric_odeint/tutorial/all_examples.html +++ /dev/null @@ -1,376 +0,0 @@ - - - -All examples - - - - - - - - -
-PrevUpHomeNext -
-
-

-All examples -

-

- The following table gives an overview over all examples. -

-
-

Table 1.5. Examples Overview

-
---- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- File -

- 		-

- Brief Description -

-
-

- harmonic_oscillator.cpp -

- 		-

- The harmonic oscillator examples gives a brief introduction to - odeint and shows the usage of the classical Runge-Kutta-solvers. -

-
-

- simple1d.cpp -

- 		-

- Integrating a simple, one-dimensional ODE showing the usage of - integrate- and generate-functions. -

-
-

- elliptic_functions.cpp -

- 		-

- Example calculating the elliptic functions using Bulirsch-Stoer - and Runge-Kutta-Dopri5 Steppers with dense output. -

-
-

- solar_system.cpp -

- 		-

- The solar system example shows the usage of the symplectic solvers. -

-
-

- chaotic_system.cpp -

- 		-

- The chaotic system examples integrates the Lorenz system and calculates - the Lyapunov exponents. -

-
-

- stiff_system.cpp -

- 		-

- The stiff system example shows the usage of the stiff solvers using - the Jacobian of the system function. -

-
-

- van_der_pol_stiff.cpp -

- 		-

- This stiff system example again shows the usage of the stiff solvers - by integrating the van der Pol osscillator. -

-
-

- stuart_landau.cpp -

- 		-

- The Stuart-Landau example shows how odeint can be used with complex - state types. -

-
-

- fpu.cpp -

- 		-

- The Fermi-Pasta-Ulam (FPU) example shows how odeint can be used - to integrate lattice systems. -

-
-

- phase_oscillator_ensemble.cpp -

- 		-

- The phase oscillator ensemble example shows how globally coupled - oscillators can be analyzed and how statistical measures can be - computed during integration. -

-
-

- stepper_details.cpp -

- 		-

- Trivial example showing the usability of the several stepper classes. -

-
-

- harmonic_oscillator_units.cpp -

- 		-

- This examples shows how Boost.Units - can be used with odeint. -

-
-

- two_dimensional_phase_lattice.cpp -

- 		-

- The 2D phase oscillator example shows how a two-dimensional lattice - works with odeint and how matrix types can be used as state types - in odeint. -

-
-

- resizing_lattice.cpp -

- 		-

- Shows the strength of odeint's memory management by simulating - a Hamiltonian system on an expanding lattice. -

-
-

- list_lattice.cpp -

- 		-

- Example of a phase lattice integration using std::list - as state type. -

-
-

- lorenz_point.cpp -

- 		-

- Alternative way of integrating lorenz by using a self defined point3d - data type as state type. -

-
-

- lorenz_gmpxx.cpp -

- 		-

- This examples integrates the Lorenz system by means of an arbitrary - precision type. -

-
-

- my_vector.cpp -

- 		-

- Simple example showing how to get odeint to work with a self-defined - vector type. -

-
-

- phase_oscillator_ensemble.cu -

- 		-

- The Thrust phase oscillator ensemble example shows how globally - coupled oscillators can be analyzed with Thrust and CUDA, employing - the power of modern graphic devices. -

-
-

- phase_oscillator_chain.cu -

- 		-

- The Thrust phase oscillator chain example shows how chains of nearest - neighbor coupled oscillators can be integrated with Thrust and - odeint. -

-
-

- lorenz_parameters.cu -

- 		-

- The Lorenz parameters examples show how ensembles of ordinary differential - equations can be solved by means of Thrust to study the dependence - of an ODE on some parameters. -

-
-

- gauss_packet.cpp -

- 		-

- The MTL-Gauss-packet example shows how the MTL can be easily used - with odeint. -

-
-

- stochastic_euler.cpp -

- 		-

- Implementation of a custom stepper - the stochastic euler - for - solving stochastic differential equations. -

-
-

- generation_functions.cpp -

- 		-

- Shows skeletal code on how to implemente own factory functions. -

-
-
-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html b/doc/boost_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html deleted file mode 100644 index 65f42787..00000000 --- a/doc/boost_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html +++ /dev/null @@ -1,278 +0,0 @@ - - - -Chaotic systems and Lyapunov exponents - - - - - - - - -
-PrevUpHomeNext -
-
-

-Chaotic - systems and Lyapunov exponents -

-

- odeint can easily be used to investigate the properties of chaotic deterministic - systems. In mathematical terms chaotic refers to an exponential growth of - perturbations δ x. In order to observe this exponential - growth one usually solves the equations for the tangential dynamics which - is again an ordinary differential equation. These equations are linear but - time dependent and can be obtained via -

-

- d δ x / dt = J(x) δ x -

-

- where J is the Jacobian of the system under consideration. - δ x can also be interpreted as a perturbation of the original - system. In principle n of these perturbations exit, - they form a hypercube and evolve in the time. The Lyapunov exponents are - then defined as logarithmic growth rates of the perturbations. If one Lyapunov - exponent is larger then zero nearby trajectories diverge exponentially hence - they are chaotic. If the largest Lyapunov exponent is zero one is usually - faced with periodic motion. In the case of a largest Lyapunov exponent smaller - then zero the converges the a fixed point. More information's about Lyapunov - exponents and nonlinear dynamical systems can be found in many textbooks, - see for example: E. Ott "Chaos is Dynamical Systems", Cambridge. -

-

- To calculate the Lyapunov exponents numerically one usually solves the equations - of motion for n perturbations and orthonormalizes them - every k steps. The Lyapunov exponent is the average - the logarithm of the stretching factor of each perturbation. -

-

- To demonstrate how one can use odeint to determine the Lyapunov exponents - we choose the Lorenz system. It is one of the most studied dynamical systems - in the nonlinear dynamics community. For the standard parameters it possesses - a strange attractor with non-integer dimension. The Lyapunov exponents take - values of approximately 0.9, 0 and -12. -

-

- The implementation of the Lorenz system is -

-

-

-
const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-typedef boost::array< double , 3 > lorenz_state_type;
-
-void lorenz( const lorenz_state_type &x , lorenz_state_type &dxdt , double t )
-{
-    dxdt[0] = sigma * ( x[1] - x[0] );
-    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
-    dxdt[2] = -b * x[2] + x[0] * x[1];
-}
-
-

- We need also to integrate the set of the perturbations. This is done in parallel - to the original system, hence within one system function. Of course, we want - to use the above definition of the Lorenz system, hence the definition of - the system function including the Lorenz system itself and the perturbation - could look like: -

-

-

-
const size_t n = 3;
-const size_t num_of_lyap = 3;
-const size_t N = n + n*num_of_lyap;
-
-typedef std::tr1::array< double , N > state_type;
-typedef std::tr1::array< double , num_of_lyap > lyap_type;
-
-void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
-{
-    lorenz( x , dxdt , t );
-
-    for( size_t l=0 ; l<num_of_lyap ; ++l )
-    {
-        const double *pert = x.begin() + 3 + l * 3;
-        double *dpert = dxdt.begin() + 3 + l * 3;
-        dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
-        dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
-        dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
-    }
-}
-
-

-

-

- The perturbations are stored linearly in the state_type - behind the state of the Lorenz system. The problem that lorenz() and lorenz_with_lyap() have different - state types. A simple trick to over come this problem is put the Lorenz system - inside a functor with templatized arguments: -

-

-

-
struct lorenz
-{
-    template< class StateIn , class StateOut , class Value >
-    void operator()( const StateIn &x , StateOut &dxdt , Value t )
-    {
-        dxdt[0] = sigma * ( x[1] - x[0] );
-        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
-        dxdt[2] = -b * x[2] + x[0] * x[1];
-    }
-};
-
-void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
-{
-    lorenz()( x , dxdt , t );
-    ...
-}
-
-

- This works fine and lorenz_with_lyap - can be used for example via -

-
state_type x;
-// initialize x
-explicit_rk4< state_type > rk4;
-integrate_n_steps( rk4 , lorenz_with_lyap , x , 0.0 , 0.01 , 1000 );
-
-

- This code snippet performs 1000 steps with constant step size 0.01. -

-

- A real world use case for the calculation of the Lyapunov exponents of Lorenz - system would always include some transient steps, just to ensure that the - current state lies on the attractor, hence it would look like -

-

-

-
state_type x;
-// initialize x
-explicit_rk4< state_type > rk4;
-integrate_n_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 );
-
-

- The problem is now, that x - is the full state containing also the perturbations and integrate_n_steps - does not know that it should only use 3 elements. In detail, odeint and its - steppers determine the length of the system under consideration by determining - the length of the state. In the classical solvers, e.g. from Numerical Recipes, - the problem was solved by pointer to the state and an appropriate length, - something similar to -

-

-

-
void lorenz( double* x , double *dxdt , double t, void* params )
-{
-    ...
-}
-
-int system_length = 3;
-rk4( x , system_length , t , dt , lorenz );
-
-

-

-

- But odeint supports a similar and much more sophisticated concept: Boost.Range. - To make the steppers and the system ready to work with Boost.Range - the system has to by changed: -

-

-

-
struct lorenz
-{
-    template< class State , class Deriv >
-    void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
-    {
-        typename boost::range_iterator< const State >::type x = boost::begin( x_ );
-        typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
-
-        dxdt[0] = sigma * ( x[1] - x[0] );
-        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
-        dxdt[2] = -b * x[2] + x[0] * x[1];
-    }
-};
-
-

-

-

- This is in principle all. Now, we only have to call integrate_n_steps - with a range including only the first 3 components of x: -

-

-

-
// perform 10000 transient steps
-integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );
-
-

-

-

- Having integrated a sufficient number of transients steps we are now able - to calculate the Lyapunov exponents: -

-
    -
  1. - We initialize the perturbations. They are stored linearly behind the - state of the Lorenz system. The perturbations are initialized such that - p ​ij = δ ​ij, where p ​ij is the - j-component of the i.-th perturbation - and δ ​ij is the Kronecker symbol. -
    2. -
  2. - Integrate 100 steps of the full system with perturbations -
    3. -
  3. - Orthonormalize the perturbation using Gram-Schmidt orthonormalization - algorithm. -
    4. -
  4. - Repeat step 2 and 3. Every 10000 steps write the current Lyapunov exponent. -
    5. -
-

-

-
fill( x.begin()+n , x.end() , 0.0 );
-for( size_t i=0 ; i<num_of_lyap ; ++i ) x[n+n*i+i] = 1.0;
-fill( lyap.begin() , lyap.end() , 0.0 );
-
-double t = 0.0;
-size_t count = 0;
-while( true )
-{
-
-    t = integrate_n_steps( rk4 , lorenz_with_lyap , x , t , dt , 100 );
-    gram_schmidt< num_of_lyap >( x , lyap , n );
-    ++count;
-
-    if( !(count % 100000) )
-    {
-        cout << t;
-        for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
-        cout << endl;
-    }
-}
-
-

-

-

- The full code can be found here: chaotic_system.cpp -

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/harmonic_oscillator.html b/doc/boost_numeric_odeint/tutorial/harmonic_oscillator.html deleted file mode 100644 index 16539c58..00000000 --- a/doc/boost_numeric_odeint/tutorial/harmonic_oscillator.html +++ /dev/null @@ -1,1489 +0,0 @@ - - - -Harmonic oscillator - - - - - - - - -
-PrevUpHomeNext -
-
-

-Harmonic - oscillator -

-
-
Define - the ODE
-
Stepper - Types
-
Integration - with Constant Step Size
-
Integration - with Adaptive Step Size
-
Using - iterators
-
-
-

-Define - the ODE -

-

- First of all, you have to specify the data type that represents a state - x of your system. Mathematically, this usually is - an n-dimensional vector with real numbers or complex numbers as scalar - objects. For odeint the most natural way is to use vector< double > or vector< complex< double > > - to represent the system state. However, odeint can deal with other container - types as well, e.g. boost::array< double , N >, as long as it is fulfills some requirements - defined below. -

-

- To integrate a differential equation numerically, one also has to define - the rhs of the equation x' = f(x). In odeint you supply - this function in terms of an object that implements the ()-operator with - a certain parameter structure. Hence, the straight forward way would be - to just define a function, e.g: -

-

-

-
/* The type of container used to hold the state vector */
-typedef std::vector< double > state_type;
-
-const double gam = 0.15;
-
-/* The rhs of x' = f(x) */
-void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
-{
-    dxdt[0] = x[1];
-    dxdt[1] = -x[0] - gam*x[1];
-}
-
-

-

-

- The parameters of the function must follow the example above where x is the current state, here a two-component - vector containing position q and momentum p - of the oscillator, dxdt - is the derivative x' and should be filled by the function - with f(x), and t - is the current time. Note that in this example t is - not required to calculate f, however odeint expects - the function signature to have exactly three parameters (there are exception, - discussed later). -

-

- A more sophisticated approach is to implement the system as a class where - the rhs function is defined as the ()-operator of the class with the same - parameter structure as above: -

-

-

-
/* The rhs of x' = f(x) defined as a class */
-class harm_osc {
-
-    double m_gam;
-
-public:
-    harm_osc( double gam ) : m_gam(gam) { }
-
-    void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
-    {
-        dxdt[0] = x[1];
-        dxdt[1] = -x[0] - m_gam*x[1];
-    }
-};
-
-

-

-

- odeint can deal with instances of such classes instead of pure functions - which allows for cleaner code. -

-
-
-

-Stepper - Types -

-

- Numerical integration works iteratively, that means you start at a state - x(t) and perform a time-step of length dt - to obtain the approximate state x(t+dt). There exist - many different methods to perform such a time-step each of which has a - certain order q. If the order of a method is q - than it is accurate up to term ~dtq that means the - error in x made by such a step is ~dtq+1. - odeint provides several steppers of different orders from which you can - choose: -

-
-

Table 1.2. Stepper Algorithms

-
----------- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Algorithm -

- 		-

- Class -

- 		-

- Concept -

- 		-

- System Concept -

- 		-

- Order -

- 		-

- Error Estimation -

- 		-

- Dense Output -

- 		-

- Internal state -

- 		-

- Remarks -

-
-

- Explicit Euler -

- 		-

- euler -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- 1 -

- 		-

- No -

- 		-

- Yes -

- 		-

- No -

- 		-

- Very simple, only for demonstrating purpose -

-
-

- Modified Midpoint -

- 		-

- modified_midpoint -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable (2) -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Used in Bulirsch-Stoer implementation -

-
-

- Runge-Kutta 4 -

- 		-

- runge_kutta4 -

- 		-

- Stepper -

- 		-

- System -

- 		-

- 4 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- The classical Runge Kutta scheme, good general scheme without - error control -

-
-

- Cash-Karp -

- 		-

- runge_kutta_cash_karp54 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 5 -

- 		-

- Yes (4) -

- 		-

- No -

- 		-

- No -

- 		-

- Good general scheme with error estimation, to be used in controlled_error_stepper -

-
-

- Dormand-Prince 5 -

- 		-

- runge_kutta_dopri5 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 5 -

- 		-

- Yes (4) -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- Standard method with error control and dense output, to be used - in controlled_error_stepper and in dense_output_controlled_explicit_fsal. -

-
-

- Fehlberg 78 -

- 		-

- runge_kutta_fehlberg78 -

- 		-

- Error - Stepper -

- 		-

- System -

- 		-

- 8 -

- 		-

- Yes (7) -

- 		-

- No -

- 		-

- No -

- 		-

- Good high order method with error estimation, to be used in controlled_error_stepper. -

-
-

- Adams Bashforth -

- 		-

- adams_bashforth -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Multistep method -

-
-

- Adams Moulton -

- 		-

- adams_moulton -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Multistep method -

-
-

- Adams Bashforth Moulton -

- 		-

- adams_bashforth_moulton -

- 		-

- Stepper -

- 		-

- System -

- 		-

- configurable -

- 		-

- No -

- 		-

- No -

- 		-

- Yes -

- 		-

- Combined multistep method -

-
-

- Controlled Runge Kutta -

- 		-

- controlled_runge_kutta -

- 		-

- Controlled - Stepper -

- 		-

- System -

- 		-

- depends -

- 		-

- Yes -

- 		-

- No -

- 		-

- depends -

- 		-

- Error control for Error - Stepper. Requires an Error - Stepper from above. Order depends on the given ErrorStepper -

-
-

- Dense Output Runge Kutta -

- 		-

- dense_output_runge_kutta -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- depends -

- 		-

- No -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- Dense output for Stepper - and Error - Stepper from above if they provide dense output functionality - (like euler and - runge_kutta_dopri5). - Order depends on the given stepper. -

-
-

- Bulirsch-Stoer -

- 		-

- bulirsch_stoer -

- 		-

- Controlled - Stepper -

- 		-

- System -

- 		-

- variable -

- 		-

- Yes -

- 		-

- No -

- 		-

- No -

- 		-

- Stepper with step size and order control. Very good if high precision - is required. -

-
-

- Bulirsch-Stoer Dense Output -

- 		-

- bulirsch_stoer_dense_out -

- 		-

- Dense - Output Stepper -

- 		-

- System -

- 		-

- variable -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Stepper with step size and order control as well as dense output. - Very good if high precision and dense output is required. -

-
-

- Implicit Euler -

- 		-

- implicit_euler -

- 		-

- Stepper -

- 		-

- Implicit - System -

- 		-

- 1 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Basic implicit routine. Requires the Jacobian. Works only with - Boost.UBlas - vectors as state types. -

-
-

- Rosenbrock 4 -

- 		-

- rosenbrock4 -

- 		-

- Error - Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Good for stiff systems. Works only with Boost.UBlas - vectors as state types. -

-
-

- Controlled Rosenbrock 4 -

- 		-

- rosenbrock4_controller -

- 		-

- Controlled - Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Rosenbrock 4 with error control. Works only with Boost.UBlas - vectors as state types. -

-
-

- Dense Output Rosenbrock 4 -

- 		-

- rosenbrock4_dense_output -

- 		-

- Dense - Output Stepper -

- 		-

- Implicit - System -

- 		-

- 4 -

- 		-

- Yes -

- 		-

- Yes -

- 		-

- No -

- 		-

- Controlled Rosenbrock 4 with dense output. Works only with Boost.UBlas - vectors as state types. -

-
-

- Symplectic Euler -

- 		-

- symplectic_euler -

- 		-

- Stepper -

- 		-

- Symplectic - System Simple - Symplectic System -

- 		-

- 1 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Basic symplectic solver for separable Hamiltonian system -

-
-

- Symplectic RKN McLachlan -

- 		-

- symplectic_rkn_sb3a_mclachlan -

- 		-

- Stepper -

- 		-

- Symplectic - System Simple - Symplectic System -

- 		-

- 6 -

- 		-

- No -

- 		-

- No -

- 		-

- No -

- 		-

- Symplectic solver for separable Hamiltonian system with order - 6 -

-
-
-

- Some of steppers in the table above are special: Some need the Jacobian - of the ODE, others are constructed for special ODE-systems like Hamiltonian - systems. We will show typical examples and use-cases in this tutorial and - which kind of steppers should be applied. -

-
-
-

-Integration - with Constant Step Size -

-

- The basic stepper just performs one time-step and doesn't give you any - information about the error that was made (except that you know it is of - order q+1). Such steppers are used with constant step - size that should be chosen small enough to have reasonable small errors. - However, you should apply some sort of validity check of your results (like - observing conserved quantities) because you have no other control of the - error. The following example defines a basic stepper based on the classical - Runge-Kutta scheme of 4th order. The declaration of the stepper requires - the state type as template parameter. The integration can now be done by - using the integrate_const( Stepper, System, state, start_time, end_time, step_size - ) function from odeint: -

-

-

-
runge_kutta4< state_type > stepper;
-integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
-
-

-

-

- This call integrates the system defined by harmonic_oscillator - using the RK4 method from t=0 to 10 - with a step-size dt=0.01 and the initial condition - given in x. The result, - x(t=10) is stored in x - (in-place). Each stepper defines a do_step - method which can used directly. So, you write down the above example as -

-

-

-
const double dt = 0.01;
-for( double t=0.0 ; t<10.0 ; t+= dt )
-    stepper.do_step( harmonic_oscillator , x , t , dt );
-
-

-

-
-
-

-Integration - with Adaptive Step Size -

-

- To improve the numerical results and additionally minimize the computational - effort, the application of a step size control is advisable. Step size - control is realized via stepper algorithms that additionally provide an - error estimation of the applied step. odeint provides a number of such - ErrorSteppers and we will show their usage - on the example of explicit_error_rk54_ck - -- a 5th order Runge-Kutta method with 4th order error estimation and coefficients - introduced by Cash and Karp. -

-

-

-
typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
-
-

-

-

- Given the error stepper, one still needs an instance that checks the error - and adjusts the step size accordingly. In odeint, this is done by ControlledSteppers. For the runge_kutta_cash_karp54 - stepper a controlled_runge_kutta - stepper exists which can be used via -

-

-

-
typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
-controlled_stepper_type controlled_stepper;
-integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
-
-

-

-

- As above, this integrates the system defined by harmonic_oscillator, - but now using an adaptive step size method based on the Runge-Kutta Cash-Karp - 54 scheme from t=0 to 10 with - an initial step size of dt=0.01 (will be adjusted) - and the initial condition given in x. The result, x(t=10), - will also be stored in x (in-place). -

-

- In the above example an error stepper is nested in a controlled stepper. - This is a nice technique; however one drawback is that one always needs - to define both steppers. Of course, one could also write the instantiation - of the controlled stepper into the call of the integrate function but a - complete knowledge of the underlying stepper types is still necessary. - Another point is, that the error tolerances for the step size control are - not easily included into the controlled stepper. Both issues can be solved - by using make_controlled: -

-

-

-
integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) ,
-                    harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
-
-

-

-

- make_controlled can be - used with many of the steppers of odeint. The first parameter is the absolute - error tolerance eps_abs and the second is the relative - error tolerance eps_rel which is used during the integration. - The template parameter determines from which error stepper a controlled - stepper should be instantiated. An alternative syntax of make_controlled is -

-

-

-
integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) ,
-                    harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
-
-

-

-

- For the Runge Kutta controller the error made during one step is compared - with eps_abs + eps_rel * ( ax * |x| + adxdt * dt * |dxdt| ). - If the error is smaller than this value the current step is accept otherwise - it is rejected and the step size is decreased. Note, that the step size - is also increased if the error gets too small compared to the rhs of the - above relation. The full instantiation of the controlled_runge_kutta - with all parameters is therefore -

-

-

-
double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
-controlled_stepper_type controlled_stepper(
-    default_error_checker< double >( abs_err , rel_err , a_x , a_dxdt ) );
-integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
-
-

-

-

- When using make_controlled - the parameter ax and adxdt are - used with their standard values of 1. -

-

- In the tables below, one can find all steppers which are working with - make_controlled and make_dense_output which is the analog - for the dense output steppers. -

-
-

Table 1.3. Generation functions make_controlled( abs_error , rel_error , stepper - )

-
----- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-

- Stepper -

- 		-

- Result of make_controlled -

- 		-

- Remarks -

-
-

- runge_kutta_cash_karp54 -

- 		-

- controlled_runge_kutta< runge_kutta_cash_karp54 - , default_error_checker<...> > -

- 		-

- ax=1, adxdt=1 -

-
-

- runge_kutta_fehlberg78 -

- 		-

- controlled_runge_kutta< runge_kutta_fehlberg78 - , default_error_checker<...> > -

- 		-

- ax=1, adxdt=1 -

-
-

- runge_kutta_dopri5 -

- 		-

- controlled_runge_kutta< runge_kutta_dopri5 - , default_error_checker<...> > -

- 		-

- a x=1, adxdt=1 -

-
-

- rosenbrock4 -

- 		-

- rosenbrock4_controlled< rosenbrock4 - > -

- 		-

- - -

-
-
-
-

Table 1.4. Generation functions make_dense_output( abs_error , rel_error , - stepper )

-
----- - - - - - - - - - - - - - - - - - -
-

- Stepper -

- 		-

- Result of make_dense_output -

- 		-

- Remarks -

-
-

- runge_kutta_dopri5 -

- 		-

- dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5 - , default_error_checker<...> > - > -

- 		-

- a x=1, adxdt=1 -

-
-

- rosenbrock4 -

- 		-

- rosenbrock4_dense_output< rosenbrock4_controller< rosenbrock4 - > > -

- 		-

- - -

-
-
-

- When using make_controlled - or make_dense_output one - should be aware which exact type is used and how the step size control - works. -

-
-
-

-Using - iterators -

-

- odeint supports iterators for solving ODEs. That is you instantiate a pair - of iterators and instead of using the integrate routines with an appropriate - observer you put the iterators in one of the algorithm from the C++ standard - library or from Boost.Range. An example is -

-

-

-
std::for_each( make_adaptive_iterator_begin( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 ) ,
-               make_adaptive_iterator_end( controlled_stepper , harmonic_oscillator , x ) ,
-               write_state() );
-
-

-

-
-

- The full cpp file for this example can be found here: harmonic_oscillator.cpp -

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/references.html b/doc/boost_numeric_odeint/tutorial/references.html deleted file mode 100644 index 13fd9573..00000000 --- a/doc/boost_numeric_odeint/tutorial/references.html +++ /dev/null @@ -1,87 +0,0 @@ - - - -References - - - - - - - - -
-PrevUpHomeNext -
-
-

-References -

-

- General information about numerical integration of - ordinary differential equations: -

-

- [1] Press William H et al., Numerical Recipes 3rd Edition: The Art of Scientific - Computing, 3rd ed. (Cambridge University Press, 2007). -

-

- [2] Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner, Solving Ordinary - Differential Equations I: Nonstiff Problems, 2nd ed. (Springer, Berlin, 2009). -

-

- [3] Ernst Hairer and Gerhard Wanner, Solving Ordinary Differential Equations - II: Stiff and Differential-Algebraic Problems, 2nd ed. (Springer, Berlin, - 2010). -

-

- Symplectic integration of numerical integration: -

-

- [4] Ernst Hairer, Gerhard Wanner, and Christian Lubich, Geometric Numerical - Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, - 2nd ed. (Springer-Verlag Gmbh, 2006). -

-

- [5] Leimkuhler Benedict and Reich Sebastian, Simulating Hamiltonian Dynamics - (Cambridge University Press, 2005). -

-

- Special symplectic methods: -

-

- [6] Haruo Yoshida, “Construction of higher order symplectic integrators,” - Physics Letters A 150, no. 5 (November 12, 1990): 262-268. -

-

- [7] Robert I. McLachlan, “On the numerical integration of ordinary differential - equations by symmetric composition methods,” SIAM J. Sci. Comput. 16, no. - 1 (1995): 151-168. -

-

- Special systems: -

-

- [8] Fermi-Pasta-Ulam - nonlinear lattice oscillations -

-

- [9] Arkady Pikovsky, Michael Rosemblum, and Jürgen Kurths, Synchronization: - A Universal Concept in Nonlinear Sciences. (Cambridge University Press, 2001). -

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/solar_system.html b/doc/boost_numeric_odeint/tutorial/solar_system.html deleted file mode 100644 index d2b567dc..00000000 --- a/doc/boost_numeric_odeint/tutorial/solar_system.html +++ /dev/null @@ -1,300 +0,0 @@ - - - -Solar system - - - - - - - - -
-PrevUpHomeNext -
-
-

-Solar system -

-
-
Gravitation - and energy conservation
-
Define - the system function
-
-
-

-Gravitation - and energy conservation -

-

- The next example in this tutorial is a simulation of the outer solar system, - consisting of the sun, Jupiter, Saturn, Uranus, Neptune and Pluto. -

-

- solar_system -

-

- Each planet and of course the sun will be represented by mass points. The - interaction force between each object is the gravitational force which - can be written as -

-

- F​ij = -γ m​i m​j ( q​i - q​j ) / | q​i - q​j | 3 -

-

- where γ is the gravitational constant, m​i - and m​j are the masses and q​i - and q​j are the locations of the two objects. The equations - of motion are then -

-

- dq​i / dt = p​i -

-

- dp​i / dt = 1 / m​i Σ​ji F​ij -

-

- where p​i is the momenta of object i. - The equations of motion can also be derived from the Hamiltonian -

-

- H = Σ​i p​i2 / ( 2 m​i ) + Σ​j V( q​i , q​j ) -

-

- with the interaction potential V(q​i,q​j). The Hamiltonian - equations give the equations of motion -

-

- dq​i / dt = dH / dp​i -

-

- dp​i / dt = -dH / dq​i -

-

- In time independent Hamiltonian system the energy and the phase space volume - are conserved and special integration methods have to be applied in order - to ensure these conservation laws. The odeint library provides classes - for separable Hamiltonian systems, which can be written in the form H - = Σ -p​i2 / (2m​i) + H​q(q), where H​q(q) only - depends on the coordinates. Although this functional form might look a - bit arbitrary, it covers nearly all classical mechanical systems with inertia - and without dissipation, or where the equations of motion can be written - in the form dq​i / dt = p​i / m​i , dp​i / dt = - f( q​i ). -

-
- - - - - -
[Note]Note

- A short physical note: While the two-body-problem is known to be integrable, - that means it can be solved with purely analytic techniques, already - the three-body-problem is not solveable. This was found in the end of - the 19th century by H. Poincare which lead to the whole new subject of - Chaos Theory. -

-
-
-

-Define - the system function -

-

- To implement this system we define a 3D point type which will represent - the space as well as the velocity. Therefore, we use the operators from - Boost.Operators: -

-

-

-
/*the point type */
-template< class T , size_t Dim >
-class point :
-    boost::additive1< point< T , Dim > ,
-    boost::additive2< point< T , Dim  > , T ,
-    boost::multiplicative2< point< T , Dim > , T
-    > > >
-    {
-    public:
-
-        const static size_t dim = Dim;
-        typedef T value_type;
-        typedef point< value_type , dim > point_type;
-
-        // ...
-        // constructors
-
-        // ...
-        // operators
-
-    private:
-
-        T m_val[dim];
-    };
-
-    //...
-    // more operators
-
-

-

-

- The next step is to define a container type storing the values of q - and p and to define system functions. As container - type we use boost::array -

-

-

-
// we simulate 5 planets and the sun
-const size_t n = 6;
-
-typedef point< double , 3 > point_type;
-typedef boost::array< point_type , n > container_type;
-typedef boost::array< double , n > mass_type;
-
-

-

-

- The container_type is different - from the state type of the ODE. The state type of the ode is simply a - pair< - container_type , - container_type > - since it needs the information about the coordinates and the momenta. -

-

- Next we define the system's equations. As we will use a stepper that accounts - for the Hamiltonian (energy-preserving) character of the system, we have - to define the rhs different from the usual case where it is just a single - function. The stepper will make use of the separable character, which means - the system will be defined by two objects representing f(p) = - -dH/dq and g(q) = dH/dp: -

-

-

-
const double gravitational_constant = 2.95912208286e-4;
-
-struct solar_system_coor
-{
-    const mass_type &m_masses;
-
-    solar_system_coor( const mass_type &masses ) : m_masses( masses ) { }
-
-    void operator()( const container_type &p , container_type &dqdt ) const
-    {
-        for( size_t i=0 ; i<n ; ++i )
-            dqdt[i] = p[i] / m_masses[i];
-    }
-};
-
-

-

-

-

-
struct solar_system_momentum
-{
-    const mass_type &m_masses;
-
-    solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }
-
-    void operator()( const container_type &q , container_type &dpdt ) const
-    {
-        const size_t n = q.size();
-        for( size_t i=0 ; i<n ; ++i )
-        {
-            dpdt[i] = 0.0;
-            for( size_t j=0 ; j<i ; ++j )
-            {
-                point_type diff = q[j] - q[i];
-                double d = abs( diff );
-                diff *= ( gravitational_constant * m_masses[i] * m_masses[j] / d / d / d );
-                dpdt[i] += diff;
-                dpdt[j] -= diff;
-
-            }
-        }
-    }
-};
-
-

-

-

- In general a three body-system is chaotic, hence we can not expect that - arbitrary initial conditions of the system will lead to a dynamic which - is comparable with the solar system. That is we have to define proper initial - conditions, which are taken from the book of Hairer, Wannier, Lubich. -

-

- As mentioned above, we need to use some special integrators in order to - conserve phase space volume. There is a well known family of such integrators, - the so-called Runge-Kutta-Nystroem solvers, which we apply here in terms - of a symplectic_rkn_sb3a_mclachlan - stepper: -

-

-

-
typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
-const double dt = 100.0;
-
-integrate_const(
-        stepper_type() ,
-        make_pair( solar_system_coor( masses ) , solar_system_momentum( masses ) ) ,
-        make_pair( boost::ref( q ) , boost::ref( p ) ) ,
-        0.0 , 200000.0 , dt , streaming_observer( cout ) );
-
-

-

-

- These integration routine was used to produce the above sketch of the solar - system. Note, that there are two particularities in this example. First, - the state of the symplectic stepper is not container_type - but a pair of container_type. - Hence, we must pass such a pair to the integrate function. Since, we want - to pass them as references we can simply pack them into Boost.Ref. - The second point is the observer, which is called with a state type, hence - a pair of container_type. - The reference wrapper is also passed, but this is not a problem at all: -

-

-

-
struct streaming_observer
-{
-    std::ostream& m_out;
-
-    streaming_observer( std::ostream &out ) : m_out( out ) { }
-
-    template< class State >
-    void operator()( const State &x , double t ) const
-    {
-        container_type &q = x.first;
-        m_out << t;
-        for( size_t i=0 ; i<q.size() ; ++i ) m_out << "\t" << q[i];
-        m_out << "\n";
-    }
-};
-
-

-

-

- The full example can be found here: solar_system.cpp -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/special_topics.html b/doc/boost_numeric_odeint/tutorial/special_topics.html deleted file mode 100644 index 8776ef65..00000000 --- a/doc/boost_numeric_odeint/tutorial/special_topics.html +++ /dev/null @@ -1,1016 +0,0 @@ - - - -Special topics - - - - - - - - -
-PrevUpHomeNext -
-
-

-Special - topics -

-
-
Complex - state types
-
Lattice - systems
-
Ensembles - of oscillators
-
Using - boost::units
-
Using - matrices as state types
-
Using - arbitrary precision floating point types
-
Self - expanding lattices
-
-
-

-Complex - state types -

-

- Thus far we have seen several examples defined for real values. Of course, - odeint can handle complex state types, hence ODEs which are defined on - complex vector spaces, as well. An example is the Stuart-Landau oscillator -

-

- d Ψ / dt = ( 1 + i η ) Ψ + ( 1 + i α ) | Ψ |2 Ψ -

-

- where Ψ and i is a complex variable. - Hence the state type of this system is complex< double >. The definition - of this ODE in C++ code is very simple -

-

-

-
typedef complex< double > state_type;
-
-struct stuart_landau
-{
-    double m_eta;
-    double m_alpha;
-
-    stuart_landau( double eta = 1.0 , double alpha = 1.0 )
-    : m_eta( eta ) , m_alpha( alpha ) { }
-
-    void operator()( const state_type &x , state_type &dxdt , double t ) const
-    {
-        const complex< double > I( 0.0 , 1.0 );
-        dxdt = ( 1.0 + m_eta * I ) * x - ( 1.0 + m_alpha * I ) * norm( x ) * x;
-    }
-};
-
-

-

-

- Of course, one can also use classical functions to implement the state - function. In this cast the Stuart-Landau oscillator looks like -

-

-

-
double eta = 1.0;
-double alpha = 1.0;
-
-void stuart_landau( const state_type &x , state_type &dxdt , double t )
-{
-    const complex< double > I( 0.0 , 1.0 );
-    dxdt[0] = ( 1.0 + m_eta * I ) * x[0] - ( 1.0 + m_alpha * I ) * norm( x[0] ) * x[0];
-}
-
-

-

-

- 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. -

-

- When chosing the stepper type one has to account for the "unusual" - state type: it is a single complex<double> opposed to the vector types used in - the previous examples. This means that no iterations over vector elements - have to be performed inside the stepper algorithm. You can enforce this - by supplying additional template arguments to the stepper including the - vector_space_algebra. Details - on the usage of algebras can be found in the section Adapt - you own state types. -

-

-

-
state_type x = complex< double >( 1.0 , 0.0 );
-
-const double dt = 0.1;
-
-typedef runge_kutta4< state_type , double , state_type , double ,
-                      vector_space_algebra > stepper_type;
-
-integrate_const( stepper_type() , stuart_landau( 2.0 , 1.0 ) , x , 0.0 , 10.0 , dt , streaming_observer( cout ) );
-
-

-

-

- The full cpp file for the Stuart Landau example can be found here stuart_landau.cpp -

-
- - - - - -
[Note]Note

- The fact that we have to configure a different algebra is solely due - to the fact that we use a non-vector state type and not to the usage - of complex values. So for, e.g. vector< complex<double> >, - this would not be required. -

-
-
-

-Lattice - systems -

-

- odeint can also be used to solve ordinary differential equations defined - on lattices. A prominent example is the Fermi-Pasta-Ulam system [8]. It - is a Hamiltonian system of nonlinear coupled harmonic oscillators. The - Hamiltonian is -

-

- H = Σ​i p​i2/2 + 1/2 ( q​i+1 - q​i )^2 + β / 4 ( q​i+1 - q​i )^4 -

-

- Remarkably, the Fermi-Pasta-Ulam system was the first numerical experiment - which has been implemented on a computer in 1953. It was studied at Los - Alamos on one of the first computer (a MANIAC I) and it triggered a whole - new tree of mathematical and physical science. -

-

- Like the 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​i / dt = p​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 -

-

-

-
typedef vector< double > container_type;
-
-struct fpu
-{
-    const double m_beta;
-
-    fpu( const double beta = 1.0 ) : m_beta( beta ) { }
-
-    // system function defining the ODE
-    void operator()( const container_type &q , container_type &dpdt ) const
-    {
-        size_t n = q.size();
-        double tmp = q[0] - 0.0;
-        double tmp2 = tmp + m_beta * tmp * tmp * tmp;
-        dpdt[0] = -tmp2;
-        for( size_t i=0 ; i<n-1 ; ++i )
-        {
-            tmp = q[i+1] - q[i];
-            tmp2 = tmp + m_beta * tmp * tmp * tmp;
-            dpdt[i] += tmp2;
-            dpdt[i+1] = -tmp2;
-        }
-        tmp = - q[n-1];
-        tmp2 = tmp + m_beta * tmp * tmp * tmp;
-        dpdt[n-1] += tmp2;
-    }
-
-    // calculates the energy of the system
-    double energy( const container_type &q , const container_type &p ) const
-    {
-        // ...
-    }
-
-    // calculates the local energy of the system
-    void local_energy( const container_type &q , const container_type &p , container_type &e ) const
-    {
-        // ...
-    }
-};
-
-

-

-

- Of course, you can also use boost::array< double , N > for the state type. -

-

- Now, you have to define your initial values and perform the integration. - All this can be easily done with the following piece of code: -

-

-

-
const size_t n = 64;
-container_type q( n , 0.0 ) , p( n , 0.0 );
-
-for( size_t i=0 ; i<n ; ++i )
-{
-    p[i] = 0.0;
-    q[i] = 32.0 * sin( double( i + 1 ) / double( n + 1 ) * M_PI );
-}
-
-
-const double dt = 0.1;
-
-typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
-fpu fpu_instance( 8.0 );
-
-integrate_const( stepper_type() , fpu_instance ,
-        make_pair( boost::ref( q ) , boost::ref( p ) ) ,
-        0.0 , 1000.0 , dt , streaming_observer( cout , fpu_instance , 10 ) );
-
-

-

-

- The observer uses a reference to the system object to calculate the local - energies: -

-

-

-
struct streaming_observer
-{
-    std::ostream& m_out;
-    const fpu &m_fpu;
-    size_t m_write_every;
-    size_t m_count;
-
-    streaming_observer( std::ostream &out , const fpu &f , size_t write_every = 100 )
-    : m_out( out ) , m_fpu( f ) , m_write_every( write_every ) , m_count( 0 ) { }
-
-    template< class State >
-    void operator()( const State &x , double t )
-    {
-        if( ( m_count % m_write_every ) == 0 )
-        {
-            container_type &q = x.first;
-            container_type &p = x.second;
-            container_type energy( q.size() );
-            m_fpu.local_energy( q , p , energy );
-            for( size_t i=0 ; i<q.size() ; ++i )
-            {
-                m_out << t << "\t" << i << "\t" << q[i] << "\t" << p[i] << "\t" << energy[i] << "\n";
-            }
-            m_out << "\n";
-            clog << t << "\t" << accumulate( energy.begin() , energy.end() , 0.0 ) << "\n";
-        }
-        ++m_count;
-    }
-};
-
-

-

-

- The full cpp file for this FPU example can be found here fpu.cpp -

-
-
-

-Ensembles - of oscillators -

-

- Another import high dimensional system of coupled ordinary differential - equations is an ensemble of N all-to-all coupled phase - oscillators [9]. It is defined as -

-

- dφ​k / dt = ω​k + ε / N Σ​j sin( φ​j - φ​k ) -

-

- The natural frequencies ω​i of each oscillator follow - some distribution and ε is the coupling strength. We - choose here a Lorentzian distribution for ω​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 -

-

- Z = K ei Θ = 1 / N Σ​kei φ​k -

-

- The definition of the system function is now a bit more complex since we - also need to store the individual frequencies of each oscillator. -

-

-

-
typedef vector< double > container_type;
-
-
-pair< double , double > calc_mean_field( const container_type &x )
-{
-    size_t n = x.size();
-    double cos_sum = 0.0 , sin_sum = 0.0;
-    for( size_t i=0 ; i<n ; ++i )
-    {
-        cos_sum += cos( x[i] );
-        sin_sum += sin( x[i] );
-    }
-    cos_sum /= double( n );
-    sin_sum /= double( n );
-
-    double K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum );
-    double Theta = atan2( sin_sum , cos_sum );
-
-    return make_pair( K , Theta );
-}
-
-
-struct phase_ensemble
-{
-    container_type m_omega;
-    double m_epsilon;
-
-    phase_ensemble( const size_t n , double g = 1.0 , double epsilon = 1.0 )
-    : m_omega( n , 0.0 ) , m_epsilon( epsilon )
-    {
-        create_frequencies( g );
-    }
-
-    void create_frequencies( double g )
-    {
-        boost::mt19937 rng;
-        boost::cauchy_distribution<> cauchy( 0.0 , g );
-        boost::variate_generator< boost::mt19937&, boost::cauchy_distribution<> > gen( rng , cauchy );
-        generate( m_omega.begin() , m_omega.end() , gen );
-    }
-
-    void set_epsilon( double epsilon ) { m_epsilon = epsilon; }
-
-    double get_epsilon( void ) const { return m_epsilon; }
-
-    void operator()( const container_type &x , container_type &dxdt , double /* t */ ) const
-    {
-        pair< double , double > mean = calc_mean_field( x );
-        for( size_t i=0 ; i<x.size() ; ++i )
-            dxdt[i] = m_omega[i] + m_epsilon * mean.first * sin( mean.second - x[i] );
-    }
-};
-
-

-

-

- 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 ε. -

-

-

-
struct statistics_observer
-{
-    double m_K_mean;
-    size_t m_count;
-
-    statistics_observer( void )
-    : m_K_mean( 0.0 ) , m_count( 0 ) { }
-
-    template< class State >
-    void operator()( const State &x , double t )
-    {
-        pair< double , double > mean = calc_mean_field( x );
-        m_K_mean += mean.first;
-        ++m_count;
-    }
-
-    double get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / double( m_count ) : 0.0 ; }
-
-    void reset( void ) { m_K_mean = 0.0; m_count = 0; }
-};
-
-

-

-

- Now, we do several integrations for different values of ε - 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 ε = - 2. -

-

-

-
const size_t n = 16384;
-const double dt = 0.1;
-
-container_type x( n );
-
-boost::mt19937 rng;
-boost::uniform_real<> unif( 0.0 , 2.0 * M_PI );
-boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
-
-// gamma = 1, the phase transition occurs at epsilon = 2
-phase_ensemble ensemble( n , 1.0 );
-statistics_observer obs;
-
-for( double epsilon = 0.0 ; epsilon < 5.0 ; epsilon += 0.1 )
-{
-    ensemble.set_epsilon( epsilon );
-    obs.reset();
-
-    // start with random initial conditions
-    generate( x.begin() , x.end() , gen );
-
-    // calculate some transients steps
-    integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 10.0 , dt );
-
-    // integrate and compute the statistics
-    integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 100.0 , dt , boost::ref( obs ) );
-    cout << epsilon << "\t" << obs.get_K_mean() << endl;
-}
-
-

-

-

- The full cpp file for this example can be found here phase_oscillator_ensemble.cpp -

-
-
-

-Using - boost::units -

-

- odeint also works well with Boost.Units - - a library for compile time unit and dimension analysis. Ot works by decoding - unit information into the types of values. For a one-dimensional unit you - can just use the Boost.Unit types as state type, deriv type and time type - and hand the vector_space_algebra - to the stepper definition and everything works just fine: -

-

-

-
typedef units::quantity< si::time , double > time_type;
-typedef units::quantity< si::length , double > length_type;
-typedef units::quantity< si::velocity , double > velocity_type;
-
-typedef runge_kutta4< length_type , double , velocity_type , time_type ,
-                      vector_space_algebra > stepper_type;
-
-

-

-

- If you want to solve more-dimensional problems the individual entries typically - have different units. That means that the state_type - is now possibly heterogeneous, meaning that every entry might have a different - type. To solve this problem, compile-time sequences from Boost.Fusion - can be used. -

-

- To illustrate how odeint works with Boost.Units - we use the harmonic oscillator as primary example. We start with defining - all quantities -

-

-

-
#include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/algebra/fusion_algebra.hpp>
-
-#include <boost/units/systems/si/length.hpp>
-#include <boost/units/systems/si/time.hpp>
-#include <boost/units/systems/si/velocity.hpp>
-#include <boost/units/systems/si/acceleration.hpp>
-#include <boost/units/systems/si/io.hpp>
-
-#include <boost/fusion/container.hpp>
-
-using namespace std;
-using namespace boost::numeric::odeint;
-namespace fusion = boost::fusion;
-namespace units = boost::units;
-namespace si = boost::units::si;
-
-typedef units::quantity< si::time , double > time_type;
-typedef units::quantity< si::length , double > length_type;
-typedef units::quantity< si::velocity , double > velocity_type;
-typedef units::quantity< si::acceleration , double > acceleration_type;
-typedef units::quantity< si::frequency , double > frequency_type;
-
-typedef fusion::vector< length_type , velocity_type > state_type;
-typedef fusion::vector< velocity_type , acceleration_type > deriv_type;
-
-

-

-

- Note, that the state_type - and the deriv_type are - now a compile-time fusion sequences. deriv_type - represents x' and is now different from the state - type as it has different unit definitions. Next, we define the ordinary - differential equation which is completely equivalent to the example in - Harmonic - Oscillator: -

-

-

-
struct oscillator
-{
-    frequency_type m_omega;
-
-    oscillator( const frequency_type &omega = 1.0 * si::hertz ) : m_omega( omega ) { }
-
-    void operator()( const state_type &x , deriv_type &dxdt , time_type t ) const
-    {
-        fusion::at_c< 0 >( dxdt ) = fusion::at_c< 1 >( x );
-        fusion::at_c< 1 >( dxdt ) = - m_omega * m_omega * fusion::at_c< 0 >( x );
-    }
-};
-
-

-

-

- Next, we instantiate an appropriate stepper. We must explicitly parametrize - the stepper with the state_type, - deriv_type, time_type. Furthermore, the iteration - over vector elements is now done by the fusion_algebra - which must also be given. For more on the state types / algebras see chapter - Adapt - you own state types. -

-

-

-
typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type , fusion_algebra > stepper_type;
-
-state_type x( 1.0 * si::meter , 0.0 * si::meter_per_second );
-
-integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , oscillator( 2.0 * si::hertz ) , x , 0.0 * si::second , 100.0 * si::second , 0.1 * si::second , streaming_observer( cout ) );
-
-

-

-

- It is quite easy but the compilation time might take very long. Furthermore, - the observer is defined a bit different -

-

-

-
struct streaming_observer
-{
-    std::ostream& m_out;
-
-    streaming_observer( std::ostream &out ) : m_out( out ) { }
-
-    struct write_element
-    {
-        std::ostream &m_out;
-        write_element( std::ostream &out ) : m_out( out ) { };
-
-        template< class T >
-        void operator()( const T &t ) const
-        {
-            m_out << "\t" << t;
-        }
-    };
-
-    template< class State , class Time >
-    void operator()( const State &x , const Time &t ) const
-    {
-        m_out << t;
-        fusion::for_each( x , write_element( m_out ) );
-        m_out << "\n";
-    }
-};
-
-

-

-
- - - - - -
[Caution]Caution

- Using Boost.Units - works nicely but compilation can be very time and memory consuming. For - example the unit test for Boost.Units - take up to 4 GB of memory at compilation. -

-

- The full cpp file for this example can be found here harmonic_oscillator_units.cpp. -

-
-
-

-Using - matrices as state types -

-

- 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. -

-

- 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. We like - phase oscillators, they are extremely easy and might serve for different - demonstration purposes. 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, -

-

- The definition of the system is -

-

-

-
typedef boost::numeric::ublas::matrix< double > state_type;
-
-struct two_dimensional_phase_lattice
-{
-    two_dimensional_phase_lattice( double gamma = 0.5 )
-    : m_gamma( gamma ) { }
-
-    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
-    {
-        size_t size1 = x.size1() , size2 = x.size2();
-
-        for( size_t i=1 ; i<size1-1 ; ++i )
-        {
-            for( size_t j=1 ; j<size2-1 ; ++j )
-            {
-                dxdt( i , j ) =
-                        coupling_func( x( i + 1 , j ) - x( i , j ) ) +
-                        coupling_func( x( i - 1 , j ) - x( i , j ) ) +
-                        coupling_func( x( i , j + 1 ) - x( i , j ) ) +
-                        coupling_func( x( i , j - 1 ) - x( i , j ) );
-            }
-        }
-
-        for( size_t i=0 ; i<x.size1() ; ++i ) dxdt( i , 0 ) = dxdt( i , x.size2() -1 ) = 0.0;
-        for( size_t j=0 ; j<x.size2() ; ++j ) dxdt( 0 , j ) = dxdt( 0 , x.size1() -1 ) = 0.0;
-    }
-
-    double coupling_func( double x ) const
-    {
-        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
-    }
-
-    double m_gamma;
-};
-
-

-

-

- This is in principle 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. -

-

- phase_lattice_2d_0000 phase_lattice_2d_0100 phase_lattice_2d_1000 -

-

- The full cpp for this example can be found here two_dimensional_phase_lattice.cpp. -

-
-
    -
  1. - [section Partial differential equations] -
    2. -
  2. -
    3. -
  3. - To be continued: -
    4. -
  4. -
    5. -
  5. - *Wave equation -
    6. -
  6. -
    7. -
  7. - *KdV -
    8. -
  8. -
    9. -
  9. - *Ginzburg-Landau -
    10. -
  10. -
    11. -
  11. - [endsect] -
    12. -
  12. -
    13. -
  13. - [section Ordinary differential equations on networks] -
    14. -
  14. -
    15. -
  15. - to be continued -
    16. -
  16. -
    17. -
  17. - [endsect] -
    18. -
-
-

-Using - arbitrary precision floating point types -

-

- Besides the classical floating point number like float, - double, complex< double > you can also use arbitrary precision - types, like the types from gmp - and mpfr. But you have to be - careful about instantiating any numbers. -

-

- For gmp types you have to set the default precision before any number is - instantiated. This can be easily done by calling mpf_set_default_prec( precision - ) as the first function in your - main program. Secondly, you can not use any global constant variables since - they will not be set with the default precision you have already set. -

-

- Here is a simple example: -

-

-

-
typedef mpf_class value_type;
-typedef boost::array< value_type , 3 > state_type;
-
-struct lorenz
-{
-    void operator()( const state_type &x , state_type &dxdt , value_type t ) const
-    {
-        const value_type sigma( 10.0 );
-        const value_type R( 28.0 );
-        const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
-
-        dxdt[0] = sigma * ( x[1] - x[0] );
-        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
-        dxdt[2] = -b * x[2] + x[0] * x[1];
-    }
-};
-
-

-

-

- which can be integrated: -

-

-

-
const int precision = 1024;
-mpf_set_default_prec( precision );
-
-state_type x = {{ value_type( 10.0 ) , value_type( 10.0 ) , value_type( 10.0 ) }};
-
-cout.precision( 1000 );
-integrate_const( runge_kutta4< state_type , value_type >() ,
-        lorenz() , x , value_type( 0.0 ) , value_type( 10.0 ) , value_type( value_type( 1.0 ) / value_type( 10.0 ) ) ,
-        streaming_observer( cout ) );
-
-

-

-
- - - - - -
[Caution]Caution

- The full support of arbitrary precision types depends on the functionality - they provide. For example, the types from gmp are lacking of functions - for calculating the power and arbitrary roots, hence they can not be - used with the controlled steppers. In detail, for full support the min( x , y ), - max( - x , - y ), - pow( - y , - y ) - must be callable. -

-

- The full example can be found at lorenz_gmpxx.cpp. -

-
-
-

-Self - expanding lattices -

-

- odeint supports changes of the state size during integration if a state_type - is used which can be resized, like std::vector. - The adjustment of the state's size has to be done from outside and the - stepper has to be instantiated with always_resizer - as the template argument for the resizer_type. - In this configuration, the stepper checks for changes in the state size - and adjust it's internal storage accordingly. -

-

- We exemplary show this for a Hamiltonian system of nonlinear, disordered - oscillators with nonlinear nearest neighbor coupling. -

-

- The system function is implemented in terms of a class that also provides - functions for calculating the energy. 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. -

-

-

-
typedef vector< double > coord_type;
-typedef pair< coord_type , coord_type > state_type;
-
-struct compacton_lattice
-{
-    const int m_max_N;
-    const double m_beta;
-    int m_pot_start_index;
-    vector< double > m_pot;
-
-    compacton_lattice( int max_N , double beta , int pot_start_index )
-        : m_max_N( max_N ) , m_beta( beta ) , m_pot_start_index( pot_start_index ) , m_pot( max_N )
-    {
-        srand( time( NULL ) );
-        // fill random potential with iid values from [0,1]
-        boost::mt19937 rng;
-        boost::uniform_real<> unif( 0.0 , 1.0 );
-        boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
-        generate( m_pot.begin() , m_pot.end() , gen );
-    }
-
-    void operator()( const coord_type &q , coord_type &dpdt )
-    {
-        // calculate dpdt = -dH/dq of this hamiltonian system
-        // dp_i/dt = - V_i * q_i^3 - beta*(q_i - q_{i-1})^3 + beta*(q_{i+1} - q_i)^3
-        const int N = q.size();
-        double diff = q[0] - q[N-1];
-        for( int i=0 ; i<N ; ++i )
-        {
-            dpdt[i] = - m_pot[m_pot_start_index+i] * q[i]*q[i]*q[i] -
-                    m_beta * diff*diff*diff;
-            diff = q[(i+1) % N] - q[i];
-            dpdt[i] += m_beta * diff*diff*diff;
-        }
-    }
-
-    void energy_distribution( const coord_type &q , const coord_type &p , coord_type &energies )
-    {
-        // computes the energy per lattice site normalized by total energy
-        const size_t N = q.size();
-        double en = 0.0;
-        for( size_t i=0 ; i<N ; i++ )
-        {
-            const double diff = q[(i+1) % N] - q[i];
-            energies[i] = p[i]*p[i]/2.0
-                + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i]/4.0
-                + m_beta/4.0 * diff*diff*diff*diff;
-            en += energies[i];
-        }
-        en = 1.0/en;
-        for( size_t i=0 ; i<N ; i++ )
-        {
-            energies[i] *= en;
-        }
-    }
-
-    double energy( const coord_type &q , const coord_type &p )
-    {
-        // calculates the total energy of the excitation
-        const size_t N = q.size();
-        double en = 0.0;
-        for( size_t i=0 ; i<N ; i++ )
-        {
-            const double diff = q[(i+1) % N] - q[i];
-            en += p[i]*p[i]/2.0
-                + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i] / 4.0
-                + m_beta/4.0 * diff*diff*diff*diff;
-        }
-        return en;
-    }
-
-    void change_pot_start( const int delta )
-    {
-        m_pot_start_index += delta;
-    }
-};
-
-

-

-

- The total size we allow is 1024 and we start with an initial state size - of 60. -

-

-

-
//start with 60 sites
-const int N_start = 60;
-coord_type q( N_start , 0.0 );
-q.reserve( max_N );
-coord_type p( N_start , 0.0 );
-p.reserve( max_N );
-// start with uniform momentum distribution over 20 sites
-fill( p.begin()+20 , p.end()-20 , 1.0/sqrt(20.0) );
-
-coord_type distr( N_start , 0.0 );
-distr.reserve( max_N );
-
-// create the system
-compacton_lattice lattice( max_N , beta , (max_N-N_start)/2 );
-
-//create the stepper, note that we use an always_resizer because state size might change during steps
-typedef symplectic_rkn_sb3a_mclachlan< coord_type , coord_type , double , coord_type , coord_type , double ,
-        range_algebra , default_operations , always_resizer > hamiltonian_stepper;
-hamiltonian_stepper stepper;
-hamiltonian_stepper::state_type state = make_pair( q , p );
-
-

-

-

- The lattice gets resized whenever the energy distribution comes close to - the borders distr[10] > - 1E-150, distr[distr.size()-10] - > 1E-150. - If we increase to the left, q - and p have to be rotated - because their resize function always appends at the end. Additionally, - the start index of the potential changes in this case. -

-

-

-
double t = 0.0;
-const double dt = 0.1;
-const int steps = 10000;
-for( int step = 0 ; step < steps ; ++step )
-{
-    stepper.do_step( boost::ref(lattice) , state , t , dt );
-    lattice.energy_distribution( state.first , state.second , distr );
-    if( distr[10] > 1E-150 )
-    {
-        do_resize( state.first , state.second , distr , state.first.size()+20 );
-        rotate( state.first.begin() , state.first.end()-20 , state.first.end() );
-        rotate( state.second.begin() , state.second.end()-20 , state.second.end() );
-        lattice.change_pot_start( -20 );
-        cout << t << ": resized left to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
-    }
-    if( distr[distr.size()-10] > 1E-150 )
-    {
-        do_resize( state.first , state.second , distr , state.first.size()+20 );
-        cout << t << ": resized right to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
-    }
-    t += dt;
-}
-
-

-

-

- The do_resize function - simply calls vector.resize of q - , p and distr. -

-

-

-
void do_resize( coord_type &q , coord_type &p , coord_type &distr , const int N )
-{
-    q.resize( N );
-    p.resize( N );
-    distr.resize( N );
-}
-
-

-

-

- The full example can be found in resizing_lattice.cpp -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/stiff_systems.html b/doc/boost_numeric_odeint/tutorial/stiff_systems.html deleted file mode 100644 index c3df6ab8..00000000 --- a/doc/boost_numeric_odeint/tutorial/stiff_systems.html +++ /dev/null @@ -1,163 +0,0 @@ - - - -Stiff systems - - - - - - - - -
-PrevUpHomeNext -
-
-

-Stiff systems -

-

- An important class of ordinary differential equations are so called stiff - system which are characterized by two or more time scales of different order. - Examples of such systems are found in chemical systems where reaction rates - of individual sub-reaction might differ over large ranges, for example: -

-

- d S​1 / dt = - 101 S​2 - 100 S​1 -

-

- d S​2 / dt = S​1 -

-

- To solve stiff systems efficiently using numerics the Jacobian -

-

- J = d f​i / d x​j -

-

- is needed. Here is the definition of the above example -

-

-

-
typedef boost::numeric::ublas::vector< double > vector_type;
-typedef boost::numeric::ublas::matrix< double > matrix_type;
-
-struct stiff_system
-{
-    void operator()( const vector_type &x , vector_type &dxdt , double /* t */ )
-    {
-        dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 * x[ 1 ];
-        dxdt[ 1 ] = x[ 0 ];
-    }
-};
-
-struct stiff_system_jacobi
-{
-    void operator()( const vector_type & /* x */ , matrix_type &J , const double & /* t */ , vector_type &dfdt )
-    {
-        J( 0 , 0 ) = -101.0;
-        J( 0 , 1 ) = -100.0;
-        J( 1 , 0 ) = 1.0;
-        J( 1 , 1 ) = 0.0;
-        dfdt[0] = 0.0;
-        dfdt[1] = 0.0;
-    }
-};
-
-

-

-

- The state type has to be a ublas::vector - and the matrix type must by a ublas::matrix - since the stiff integrator only accepts these types. However, you might want - use non-stiff intgrators on this system, too - we will do so later for demonstration. - Therefore we want to use the same function also with other state_types, realized - by templatizing the operator(): -

-

-

-
typedef boost::numeric::ublas::vector< double > vector_type;
-typedef boost::numeric::ublas::matrix< double > matrix_type;
-
-struct stiff_system
-{
-    template< class State >
-    void operator()( const State &x , State &dxdt , double t )
-    {
-        ...
-    }
-};
-
-struct stiff_system_jacobi
-{
-    template< class State , class Matrix >
-    void operator()( const State &x , Matrix &J , const double &t , State &dfdt )
-    {
-        ...
-    }
-};
-
-

-

-

- Now you can use stiff_system - in combination with std::vector or boost::array. - In the example the explicit time derivative of f(x,t) - is introduced separately in the Jacobian. If df / dt = 0 - simply fill dfdt with zeros. -

-

- A well know solver for stiff systems is the so called Rosenbrock method. - It has a step size control and dense output facilities and can be used like - all the other stepper: -

-

-

-
vector_type x( 2 , 1.0 );
-
-size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< double > >( 1.0e-6 , 1.0e-6 ) ,
-        make_pair( stiff_system() , stiff_system_jacobi() ) ,
-        x , 0.0 , 50.0 , 0.01 ,
-        cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
-
-

-

-

- During the integration 71 steps have been done. Comparing to a classical - Runge-Kutta solver this is a very good result. For example the Dormand-Prince - 5 method with step size control and dense output yields 1531 steps. -

-

-

-
vector_type x2( 2 , 1.0 );
-
-size_t num_of_steps2 = integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) ,
-        stiff_system() , x2 , 0.0 , 50.0 , 0.01 ,
-        cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
-
-

-

-

- Note, that we have used Boost.Phoenix, - a great functional programming library to create and compose the observer. -

-

- The full example can be found here: stiff_system.cpp -

-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boost_numeric_odeint/tutorial/using_cuda_and_thrust.html b/doc/boost_numeric_odeint/tutorial/using_cuda_and_thrust.html deleted file mode 100644 index 0f1cfe05..00000000 --- a/doc/boost_numeric_odeint/tutorial/using_cuda_and_thrust.html +++ /dev/null @@ -1,691 +0,0 @@ - - - -Using Cuda and Thrust - - - - - - - - -
-PrevUpHomeNext -
-
-

-Using - Cuda and Thrust -

-
-
Phase - oscillator ensemble
-
Large - oscillator chains
-
Parameter - studies
-
-

- Modern graphic cards (graphic processing units - GPUs) can be used to speed - up the performance of time consuming algorithms by means of massive parallelization. - They are designed to execute many operations in parallel. odeint can utilize - the power of GPUs by means of CUDA and Thrust, - which is a STL-like interface for the native CUDA API. -

-
- - - - - -
[Important]Important

- Thrust also supports parallelization using OpenMP. You can switch between - CUDA-parallelization and OpenMP-parallelization by a simple compiler switch. - Hence, this also provides an easy way to get basic OpenMP parallelization - into odeint. -

-

- To use odeint with CUDA a few points have to be taken into account. First - of all, the problem has to be well chosen. It makes absolutely no sense to - try to parallelize the code for a three dimensional system, it is simply - to small and not worth the effort. One single function call (kernel execution) - on the GPU is slow but you can do the operation on a huge set of data with - only one call. We have experienced that the vector size over which is parallelized - should be of the order of 106 to make full use of the - GPU. Secondly, you have to use Thrust's - algorithms and functors when implementing the rhs the ODE. This might be - tricky since it involves some kind of functional programming knowledge. -

-

- Typical applications for CUDA and odeint are large systems, like lattices - or discretizations of PDE, and parameter studies. We introduce now three - examples which show how the power of GPUs can be used in combination with - odeint. -

-
- - - - - -
[Important]Important

- The full power of CUDA is only available for really large systems where - the number of coupled ordinary differential equations is of order N=106 - or larger. For smaller systems the CPU is usually much faster. You can - also integrate an ensemble of different uncoupled ODEs in parallel as shown - in the last example. -

-
-

-Phase - oscillator ensemble -

-

- The first example is the phase oscillator ensemble from the previous section: -

-

- dφ​k / dt = ω​k + ε / N Σ​j sin( φ​j - φ​k ). -

-

- It has a phase transition at ε = 2 in the limit of infinite - numbers of oscillators N. In the case of finite N - this transition is smeared out but still clearly visible. -

-

- Thrust and Cuda are - perfectly suited for such kinds of problems where one needs a large number - of particles (oscillators). We start by defining the state type which is - a thrust::device_vector. The content of this vector - lives on the GPU. If you are not familiar with this we recommend reading - the Getting started section on the Thrust - website. -

-

-

-
//change this to float if your device does not support double computation
-typedef double value_type;
-
-//change this to host_vector< ... > of you want to run on CPU
-typedef thrust::device_vector< value_type > state_type;
-// typedef thrust::host_vector< value_type > state_type;
-
-

-

-

- Thrust follows a functional programming approach. If you want to perform - a calculation on the GPU you usually have to call a global function like - thrust::for_each, thrust::reduce, - ... with an appropriate local functor which performs the basic operation. - An example is -

-
struct add_two
-{
-    template< class T >
-    __host__ __device__
-    void operator()( T &t ) const
-    {
-        t += T( 2 );
-    }
-};
-
-// ...
-
-thrust::for_each( x.begin() , x.end() , add_two() );
-
-

- This code generically adds two to every element in the container x. -

-

- For the purpose of integrating the phase oscillator ensemble we need -

-
    -
  • - to calculate the system function, hence the r.h.s. of the ODE. -
    • -
  • - this involves computing the mean field of the oscillator example, i.e. - the values of R and θ -
    • -
-

- The mean field is calculated in a class mean_field_calculator -

-

-

-
struct mean_field_calculator
-{
-    struct sin_functor : public thrust::unary_function< value_type , value_type >
-    {
-        __host__ __device__
-        value_type operator()( value_type x) const
-        {
-            return sin( x );
-        }
-    };
-
-    struct cos_functor : public thrust::unary_function< value_type , value_type >
-    {
-        __host__ __device__
-        value_type operator()( value_type x) const
-        {
-            return cos( x );
-        }
-    };
-
-    static std::pair< value_type , value_type > get_mean( const state_type &x )
-    {
-        value_type sin_sum = thrust::reduce(
-                thrust::make_transform_iterator( x.begin() , sin_functor() ) ,
-                thrust::make_transform_iterator( x.end() , sin_functor() ) );
-        value_type cos_sum = thrust::reduce(
-                thrust::make_transform_iterator( x.begin() , cos_functor() ) ,
-                thrust::make_transform_iterator( x.end() , cos_functor() ) );
-
-        cos_sum /= value_type( x.size() );
-        sin_sum /= value_type( x.size() );
-
-        value_type K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum );
-        value_type Theta = atan2( sin_sum , cos_sum );
-
-        return std::make_pair( K , Theta );
-    }
-};
-
-

-

-

- Inside this class two member structures sin_functor - and cos_functor are defined. - They compute the sine and the cosine of a value and they are used within - a transform iterator to calculate the sum of sin(φ​k) - and cos(φ​k). The classifiers __host__ - and __device__ are CUDA - specific and define a function or operator which can be executed on the - GPU as well as on the CPU. The line -

-

-

-
value_type sin_sum = thrust::reduce(
-        thrust::make_transform_iterator( x.begin() , sin_functor() ) ,
-        thrust::make_transform_iterator( x.end() , sin_functor() ) );
-
-

-

-

- performs the calculation of this sine-sum on the GPU (or on the CPU, depending - on your thrust configuration). -

-

- The system function is defined via -

-

-

-
class phase_oscillator_ensemble
-{
-
-public:
-
-    struct sys_functor
-    {
-        value_type m_K , m_Theta , m_epsilon;
-
-        sys_functor( value_type K , value_type Theta , value_type epsilon )
-        : m_K( K ) , m_Theta( Theta ) , m_epsilon( epsilon ) { }
-
-        template< class Tuple >
-        __host__ __device__
-        void operator()( Tuple t )
-        {
-            thrust::get<2>(t) = thrust::get<1>(t) + m_epsilon * m_K * sin( m_Theta - thrust::get<0>(t) );
-        }
-    };
-
-    // ...
-
-    void operator() ( const state_type &x , state_type &dxdt , const value_type dt ) const
-    {
-        std::pair< value_type , value_type > mean_field = mean_field_calculator::get_mean( x );
-
-        thrust::for_each(
-                thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_omega.begin() , dxdt.begin() ) ),
-                thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_omega.end() , dxdt.end()) ) ,
-                sys_functor( mean_field.first , mean_field.second , m_epsilon )
-                );
-    }
-
-    // ...
-};
-
-

-

-

- This class is used within the do_step - and integrate method. It - defines a member structure sys_functor - for the r.h.s. of each individual oscillator and the operator() for the use in the steppers and integrators - of odeint. The functor computes first the mean field of φ​k - and secondly calculates the whole r.h.s. of the ODE using this mean field. - Note, how nicely thrust::tuple - and thrust::zip_iterator play together. -

-

- Now, we are ready to put everything together. All we have to do for making - odeint ready for using the GPU is to parametrize the stepper with the appropriate - thrust algebra/operations: -

-

-

-
typedef runge_kutta4< state_type , value_type , state_type , value_type , thrust_algebra , thrust_operations > stepper_type;
-
-

-

-

- Of course, you can also use a controlled or dense output stepper, e.g. -

-

-

-
typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type , thrust_algebra , thrust_operations > stepper_type;
-
-

-

-

- Then, it is straightforward to integrate the phase ensemble by creating - an instance of the rhs class and using an integration function: -

-

-

-
phase_oscillator_ensemble ensemble( N , 1.0 );
-
-

-

-

-

-
size_t steps1 = integrate_const( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , boost::ref( ensemble ) , x , 0.0 , t_transients , dt );
-
-

-

-

- We have to use boost::ref here in order to pass the rhs class - as reference and not by value. This ensures that the natural frequencies - of each oscillator are not copied when calling integrate_const. - In the full example the performance and results of the Runge-Kutta-4 and - the Dopri5 solver are compared. -

-

- The full example can be found at phase_oscillator_example.cu. -

-
-
-

-Large - oscillator chains -

-

- The next example is a large, one-dimensional chain of nearest-neighbor - coupled phase oscillators with the following equations of motion: -

-

- d φ​k / dt = ω​k + sin( φ​k+1 - φ​k ) + sin( φ​k - φ​k-1) -

-

- In principle we can use all the techniques from the previous phase oscillator - ensemble example, but we have to take special care about the coupling of - the oscillators. To efficiently implement the coupling you can use a very - elegant way employing Thrust's permutation iterator. A permutation iterator - behaves like a normal iterator on a vector but it does not iterate along - the usual order of the elements. It rather iterates along some permutation - of the elements defined by some index map. To realize the nearest neighbor - coupling we create one permutation iterator which travels one step behind - a usual iterator and another permutation iterator which travels one step - in front. The full system class is: -

-

-

-
//change this to host_vector< ... > if you want to run on CPU
-typedef thrust::device_vector< value_type > state_type;
-typedef thrust::device_vector< size_t > index_vector_type;
-//typedef thrust::host_vector< value_type > state_type;
-//typedef thrust::host_vector< size_t > index_vector_type;
-
-class phase_oscillators
-{
-
-public:
-
-    struct sys_functor
-    {
-        template< class Tuple >
-        __host__ __device__
-        void operator()( Tuple t )  // this functor works on tuples of values
-        {
-            // first, unpack the tuple into value, neighbors and omega
-            const value_type phi = thrust::get<0>(t);
-            const value_type phi_left = thrust::get<1>(t);  // left neighbor
-            const value_type phi_right = thrust::get<2>(t); // right neighbor
-            const value_type omega = thrust::get<3>(t);
-            // the dynamical equation
-            thrust::get<4>(t) = omega + sin( phi_right - phi ) + sin( phi - phi_left );
-        }
-    };
-
-    phase_oscillators( const state_type &omega )
-        : m_omega( omega ) , m_N( omega.size() ) , m_prev( omega.size() ) , m_next( omega.size() )
-    {
-        // build indices pointing to left and right neighbours
-        thrust::counting_iterator<size_t> c( 0 );
-        thrust::copy( c , c+m_N-1 , m_prev.begin()+1 );
-        m_prev[0] = 0; // m_prev = { 0 , 0 , 1 , 2 , 3 , ... , N-1 }
-
-        thrust::copy( c+1 , c+m_N , m_next.begin() );
-        m_next[m_N-1] = m_N-1; // m_next = { 1 , 2 , 3 , ... , N-1 , N-1 }
-    }
-
-    void operator() ( const state_type &x , state_type &dxdt , const value_type dt )
-    {
-        thrust::for_each(
-                thrust::make_zip_iterator(
-                        thrust::make_tuple(
-                                x.begin() ,
-                                thrust::make_permutation_iterator( x.begin() , m_prev.begin() ) ,
-                                thrust::make_permutation_iterator( x.begin() , m_next.begin() ) ,
-                                m_omega.begin() ,
-                                dxdt.begin()
-                                ) ),
-                thrust::make_zip_iterator(
-                        thrust::make_tuple(
-                                x.end() ,
-                                thrust::make_permutation_iterator( x.begin() , m_prev.end() ) ,
-                                thrust::make_permutation_iterator( x.begin() , m_next.end() ) ,
-                                m_omega.end() ,
-                                dxdt.end()) ) ,
-                sys_functor()
-                );
-    }
-
-private:
-
-    const state_type &m_omega;
-    const size_t m_N;
-    index_vector_type m_prev;
-    index_vector_type m_next;
-};
-
-

-

-

- Note, how easy you can obtain the value for the left and right neighboring - oscillator in the system functor using the permutation iterators. But, - the call of the thrust::for_each - function looks relatively complicated. Every term of the r.h.s. of the - ODE is resembled by one iterator packed in exactly the same way as it is - unpacked in the functor above. -

-

- Now we put everything together. We create random initial conditions and - decreasing frequencies such that we should get synchronization. We copy - the frequencies and the initial conditions onto the device and finally - initialize and perform the integration. As result we simply write out the - current state, hence the phase of each oscillator. -

-

-

-
// create initial conditions and omegas on host:
-vector< value_type > x_host( N );
-vector< value_type > omega_host( N );
-for( size_t i=0 ; i<N ; ++i )
-{
-    x_host[i] = 2.0 * pi * drand48();
-    omega_host[i] = ( N - i ) * epsilon; // decreasing frequencies
-}
-
-// copy to device
-state_type x = x_host;
-state_type omega = omega_host;
-
-// create stepper
-runge_kutta4< state_type , value_type , state_type , value_type , thrust_algebra , thrust_operations > stepper;
-
-// create phase oscillator system function
-phase_oscillators sys( omega );
-
-// integrate
-integrate_const( stepper , sys , x , 0.0 , 10.0 , dt );
-
-thrust::copy( x.begin() , x.end() , std::ostream_iterator< value_type >( std::cout , "\n" ) );
-std::cout << std::endl;
-
-

-

-

- The full example can be found at phase_oscillator_chain.cu. -

-
-
-

-Parameter - studies -

-

- Another important use case for Thrust - and Cuda are parameter studies of relatively small systems. Consider for - example the three-dimensional Lorenz system from the chaotic systems example - in the previous section which has three parameters. If you want to study - the behavior of this system for different parameters you usually have to - integrate the system for many parameter values. Using thrust and odeint - you can do this integration in parallel, hence you integrate a whole ensemble - of Lorenz systems where each individual realization has a different parameter - value. -

-

- In the following we will show how you can use Thrust - to integrate the above mentioned ensemble of Lorenz systems. We will vary - only the parameter β but it is straightforward to vary - other parameters or even two or all three parameters. Furthermore, we will - use the largest Lyapunov exponent to quantify the behavior of the system - (chaoticity). -

-

- We start by defining the range of the parameters we want to study. Of course, - the state_type is again a thrust::device_vector< value_type - >. -

-

-

-
vector< value_type > beta_host( N );
-const value_type beta_min = 0.0 , beta_max = 56.0;
-for( size_t i=0 ; i<N ; ++i )
-    beta_host[i] = beta_min + value_type( i ) * ( beta_max - beta_min ) / value_type( N - 1 );
-
-state_type beta = beta_host;
-
-

-

-

- The next thing we have to implement is the Lorenz system without perturbations. - Later, a system with perturbations is also implemented in order to calculate - the Lyapunov exponent. We will use an ansatz where each device function - calculates one particular realization of the Lorenz ensemble -

-

-

-
struct lorenz_system
-{
-    struct lorenz_functor
-    {
-        template< class T >
-        __host__ __device__
-        void operator()( T t ) const
-        {
-            // unpack the parameter we want to vary and the Lorenz variables
-            value_type R = thrust::get< 3 >( t );
-            value_type x = thrust::get< 0 >( t );
-            value_type y = thrust::get< 1 >( t );
-            value_type z = thrust::get< 2 >( t );
-            thrust::get< 4 >( t ) = sigma * ( y - x );
-            thrust::get< 5 >( t ) = R * x - y - x * z;
-            thrust::get< 6 >( t ) = -b * z + x * y ;
-
-        }
-    };
-
-    lorenz_system( size_t N , const state_type &beta )
-    : m_N( N ) , m_beta( beta ) { }
-
-    template< class State , class Deriv >
-    void operator()(  const State &x , Deriv &dxdt , value_type t ) const
-    {
-        thrust::for_each(
-                thrust::make_zip_iterator( thrust::make_tuple(
-                        boost::begin( x ) ,
-                        boost::begin( x ) + m_N ,
-                        boost::begin( x ) + 2 * m_N ,
-                        m_beta.begin() ,
-                        boost::begin( dxdt ) ,
-                        boost::begin( dxdt ) + m_N ,
-                        boost::begin( dxdt ) + 2 * m_N  ) ) ,
-                thrust::make_zip_iterator( thrust::make_tuple(
-                        boost::begin( x ) + m_N ,
-                        boost::begin( x ) + 2 * m_N ,
-                        boost::begin( x ) + 3 * m_N ,
-                        m_beta.begin() ,
-                        boost::begin( dxdt ) + m_N ,
-                        boost::begin( dxdt ) + 2 * m_N ,
-                        boost::begin( dxdt ) + 3 * m_N  ) ) ,
-                lorenz_functor() );
-    }
-    size_t m_N;
-    const state_type &m_beta;
-};
-
-

-

-

- As state_type a thrust::device_vector or a Boost.Range - of a device_vector is used. - The length of the state is 3N where N - is the number of systems. The system is encoded into this vector such that - all x components come first, then every y - components and finally every z components. Implementing - the device function is then a simple task, you only have to decompose the - tuple originating from the zip iterators. -

-

- Besides the system without perturbations we furthermore need to calculate - the system including linearized equations governing the time evolution - of small perturbations. Using the method from above this is straightforward, - with a small difficulty that Thrust's tuples have a maximal arity of 10. - But this is only a small problem since we can create a zip iterator packed - with zip iterators. So the top level zip iterator contains one zip iterator - for the state, one normal iterator for the parameter, and one zip iterator - for the derivative. Accessing the elements of this tuple in the system - function is then straightforward, you unpack the tuple with thrust::get<>(). - We will not show the code here, it is to large. It can be found here and - is easy to understand. -

-

- Furthermore, we need an observer which determines the norm of the perturbations, - normalizes them and averages the logarithm of the norm. The device functor - which is used within this observer is defined -

-

-

-
struct lyap_functor
-{
-    template< class T >
-    __host__ __device__
-    void operator()( T t ) const
-    {
-        value_type &dx = thrust::get< 0 >( t );
-        value_type &dy = thrust::get< 1 >( t );
-        value_type &dz = thrust::get< 2 >( t );
-        value_type norm = sqrt( dx * dx + dy * dy + dz * dz );
-        dx /= norm;
-        dy /= norm;
-        dz /= norm;
-        thrust::get< 3 >( t ) += log( norm );
-    }
-};
-
-

-

-

- Note, that this functor manipulates the state, i.e. the perturbations. -

-

- Now we complete the whole code to calculate the Lyapunov exponents. First, - we have to define a state vector. This vector contains 6N - entries, the state x,y,z and its perturbations dx,dy,dz. - We initialize them such that x=y=z=10, dx=1, - and dy=dz=0. We define a stepper type, a controlled - Runge-Kutta Dormand-Prince 5 stepper. We start with some integration to - overcome the transient behavior. For this, we do not involve the perturbation - and run the algorithm only on the state x,y,z without - any observer. Note, how Boost.Range - is used for partial integration of the state vector without perturbations - (the first half of the whole state). After the transient, the full system - with perturbations is integrated and the Lyapunov exponents are calculated - and written to stdout. -

-

-

-
state_type x( 6 * N );
-
-// initialize x,y,z
-thrust::fill( x.begin() , x.begin() + 3 * N , 10.0 );
-
-// initial dx
-thrust::fill( x.begin() + 3 * N , x.begin() + 4 * N , 1.0 );
-
-// initialize dy,dz
-thrust::fill( x.begin() + 4 * N , x.end() , 0.0 );
-
-
-// create error stepper, can be used with make_controlled or make_dense_output
-typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type , thrust_algebra , thrust_operations > stepper_type;
-
-
-lorenz_system lorenz( N , beta );
-lorenz_perturbation_system lorenz_perturbation( N , beta );
-lyap_observer obs( N , 1 );
-
-// calculate transients
-integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz , std::make_pair( x.begin() , x.begin() + 3 * N ) , 0.0 , 10.0 , dt );
-
-// calculate the Lyapunov exponents -- the main loop
-double t = 0.0;
-while( t < 10000.0 )
-{
-    integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz_perturbation , x , t , t + 1.0 , 0.1 );
-    t += 1.0;
-    obs( x , t );
-}
-
-vector< value_type > lyap( N );
-obs.fill_lyap( lyap );
-
-for( size_t i=0 ; i<N ; ++i )
-    cout << beta_host[i] << "\t" << lyap[i] << "\n";
-
-

-

-

- The full example can be found at lorenz_parameters.cu. -

-
-
- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHomeNext -
- - diff --git a/doc/boostbook.css b/doc/boostbook.css deleted file mode 100644 index d8f7bf4d..00000000 --- a/doc/boostbook.css +++ /dev/null @@ -1,618 +0,0 @@ -/*============================================================================= - Copyright (c) 2004 Joel de Guzman - http://spirit.sourceforge.net/ - - Distributed under the Boost Software License, Version 1.0. (See accompany- - ing file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -=============================================================================*/ - -/*============================================================================= - Body defaults -=============================================================================*/ - - body - { - margin: 1em; - font-family: sans-serif; - } - -/*============================================================================= - Paragraphs -=============================================================================*/ - - p - { - text-align: left; - font-size: 10pt; - line-height: 1.15; - } - -/*============================================================================= - Program listings -=============================================================================*/ - - /* Code on paragraphs */ - p tt.computeroutput - { - font-size: 9pt; - } - - pre.synopsis - { - font-size: 9pt; - margin: 1pc 4% 0pc 4%; - padding: 0.5pc 0.5pc 0.5pc 0.5pc; - } - - .programlisting, - .screen - { - font-size: 9pt; - display: block; - margin: 1pc 4% 0pc 4%; - padding: 0.5pc 0.5pc 0.5pc 0.5pc; - } - - /* Program listings in tables don't get borders */ - td .programlisting, - td .screen - { - margin: 0pc 0pc 0pc 0pc; - padding: 0pc 0pc 0pc 0pc; - } - -/*============================================================================= - Headings -=============================================================================*/ - - h1, h2, h3, h4, h5, h6 - { - text-align: left; - margin: 1em 0em 0.5em 0em; - font-weight: bold; - } - - h1 { font-size: 140%; } - h2 { -/* changed by Mario Mulansky to show logo instead of title */ -height: 127px; -background: white url(logo.jpg) left no-repeat; -text-indent: -999px; -/* end change */ -font-weight: bold; font-size: 140%; -} - h3 { font-weight: bold; font-size: 130%; } - h4 { font-weight: bold; font-size: 120%; } - h5 { font-weight: normal; font-style: italic; font-size: 110%; } - h6 { font-weight: normal; font-style: italic; font-size: 100%; } - - /* Top page titles */ - title, - h1.title, - h2.title - h3.title, - h4.title, - h5.title, - h6.title, - .refentrytitle - { - font-weight: bold; - margin-bottom: 1pc; - } - - h1.title { font-size: 140% } - h2.title { font-size: 140% } - h3.title { font-size: 130% } - h4.title { font-size: 120% } - h5.title { font-size: 110% } - h6.title { font-size: 100% } - - .section h1 - { - margin: 0em 0em 0.5em 0em; - font-size: 140%; - } - - .section h2 { font-size: 140% } - .section h3 { font-size: 130% } - .section h4 { font-size: 120% } - .section h5 { font-size: 110% } - .section h6 { font-size: 100% } - - /* Code on titles */ - h1 tt.computeroutput { font-size: 140% } - h2 tt.computeroutput { font-size: 140% } - h3 tt.computeroutput { font-size: 130% } - h4 tt.computeroutput { font-size: 130% } - h5 tt.computeroutput { font-size: 130% } - h6 tt.computeroutput { font-size: 130% } - - -/*============================================================================= - Author -=============================================================================*/ - - h3.author - { - font-size: 100% - } - -/*============================================================================= - Lists -=============================================================================*/ - - li - { - font-size: 10pt; - line-height: 1.3; - } - - /* Unordered lists */ - ul - { - text-align: left; - } - - /* Ordered lists */ - ol - { - text-align: left; - } - -/*============================================================================= - Links -=============================================================================*/ - - a - { - text-decoration: none; /* no underline */ - } - - a:hover - { - text-decoration: underline; - } - -/*============================================================================= - Spirit style navigation -=============================================================================*/ - - .spirit-nav - { - text-align: right; - } - - .spirit-nav a - { - color: white; - padding-left: 0.5em; - } - - .spirit-nav img - { - border-width: 0px; - } - -/*============================================================================= - Copyright footer -=============================================================================*/ - .copyright-footer - { - text-align: right; - font-size: 70%; - } - - .copyright-footer p - { - text-align: right; - font-size: 80%; - } - -/*============================================================================= - Table of contents -=============================================================================*/ - - .toc - { - margin: 1pc 4% 0pc 4%; - padding: 0.1pc 1pc 0.1pc 1pc; - font-size: 80%; - line-height: 1.15; - } - - .boost-toc - { - float: right; - padding: 0.5pc; - } - - /* Code on toc */ - .toc .computeroutput { font-size: 120% } - - /* No margin on nested menus */ - - .toc dl dl { margin: 0; } - -/*============================================================================= - Tables -=============================================================================*/ - - .table-title, - div.table p.title - { - margin-left: 4%; - padding-right: 0.5em; - padding-left: 0.5em; - } - - .informaltable table, - .table table - { - width: 92%; - margin-left: 4%; - margin-right: 4%; - } - - div.informaltable table, - div.table table - { - padding: 4px; - } - - /* Table Cells */ - div.informaltable table tr td, - div.table table tr td - { - padding: 0.5em; - text-align: left; - font-size: 9pt; - } - - div.informaltable table tr th, - div.table table tr th - { - padding: 0.5em 0.5em 0.5em 0.5em; - border: 1pt solid white; - font-size: 80%; - } - - table.simplelist - { - width: auto !important; - margin: 0em !important; - padding: 0em !important; - border: none !important; - } - table.simplelist td - { - margin: 0em !important; - padding: 0em !important; - text-align: left !important; - font-size: 9pt !important; - border: none !important; - } - -/*============================================================================= - Blurbs -=============================================================================*/ - - div.note, - div.tip, - div.important, - div.caution, - div.warning, - p.blurb - { - font-size: 9pt; /* A little bit smaller than the main text */ - line-height: 1.2; - display: block; - margin: 1pc 4% 0pc 4%; - padding: 0.5pc 0.5pc 0.5pc 0.5pc; - } - - p.blurb img - { - padding: 1pt; - } - -/*============================================================================= - Variable Lists -=============================================================================*/ - - div.variablelist - { - margin: 1em 0; - } - - /* Make the terms in definition lists bold */ - div.variablelist dl dt, - span.term - { - font-weight: bold; - font-size: 10pt; - } - - div.variablelist table tbody tr td - { - text-align: left; - vertical-align: top; - padding: 0em 2em 0em 0em; - font-size: 10pt; - margin: 0em 0em 0.5em 0em; - line-height: 1; - } - - div.variablelist dl dt - { - margin-bottom: 0.2em; - } - - div.variablelist dl dd - { - margin: 0em 0em 0.5em 2em; - font-size: 10pt; - } - - div.variablelist table tbody tr td p, - div.variablelist dl dd p - { - margin: 0em 0em 0.5em 0em; - line-height: 1; - } - -/*============================================================================= - Misc -=============================================================================*/ - - /* Title of books and articles in bibliographies */ - span.title - { - font-style: italic; - } - - span.underline - { - text-decoration: underline; - } - - span.strikethrough - { - text-decoration: line-through; - } - - /* Copyright, Legal Notice */ - div div.legalnotice p - { - text-align: left - } - -/*============================================================================= - Colors -=============================================================================*/ - - @media screen - { - body { - background-color: #FFFFFF; - color: #000000; - } - - /* Syntax Highlighting */ - .keyword { color: #0000AA; } - .identifier { color: #000000; } - .special { color: #707070; } - .preprocessor { color: #402080; } - .char { color: teal; } - .comment { color: #800000; } - .string { color: teal; } - .number { color: teal; } - .white_bkd { background-color: #FFFFFF; } - .dk_grey_bkd { background-color: #999999; } - - /* Links */ - a, a .keyword, a .identifier, a .special, a .preprocessor - a .char, a .comment, a .string, a .number - { - color: #005a9c; - } - - a:visited, a:visited .keyword, a:visited .identifier, - a:visited .special, a:visited .preprocessor a:visited .char, - a:visited .comment, a:visited .string, a:visited .number - { - color: #9c5a9c; - } - - h1 a, h2 a, h3 a, h4 a, h5 a, h6 a, - h1 a:hover, h2 a:hover, h3 a:hover, h4 a:hover, h5 a:hover, h6 a:hover, - h1 a:visited, h2 a:visited, h3 a:visited, h4 a:visited, h5 a:visited, h6 a:visited - { - text-decoration: none; /* no underline */ - color: #000000; - } - - /* Copyright, Legal Notice */ - .copyright - { - color: #666666; - font-size: small; - } - - div div.legalnotice p - { - color: #666666; - } - - /* Program listing */ - pre.synopsis - { - border: 1px solid #DCDCDC; - } - - .programlisting, - .screen - { - border: 1px solid #DCDCDC; - } - - td .programlisting, - td .screen - { - border: 0px solid #DCDCDC; - } - - /* Blurbs */ - div.note, - div.tip, - div.important, - div.caution, - div.warning, - p.blurb - { - border: 1px solid #DCDCDC; - } - - /* Table of contents */ - .toc - { - border: 1px solid #DCDCDC; - } - - /* Tables */ - div.informaltable table tr td, - div.table table tr td - { - border: 1px solid #DCDCDC; - } - - div.informaltable table tr th, - div.table table tr th - { - background-color: #F0F0F0; - border: 1px solid #DCDCDC; - } - - .copyright-footer - { - color: #8F8F8F; - } - - /* Misc */ - span.highlight - { - color: #00A000; - } - } - - @media print - { - /* Links */ - a - { - color: black; - } - - a:visited - { - color: black; - } - - .spirit-nav - { - display: none; - } - - /* Program listing */ - pre.synopsis - { - border: 1px solid gray; - } - - .programlisting, - .screen - { - border: 1px solid gray; - } - - td .programlisting, - td .screen - { - border: 0px solid #DCDCDC; - } - - /* Table of contents */ - .toc - { - border: 1px solid gray; - } - - .informaltable table, - .table table - { - border: 1px solid gray; - border-collapse: collapse; - } - - /* Tables */ - div.informaltable table tr td, - div.table table tr td - { - border: 1px solid gray; - } - - div.informaltable table tr th, - div.table table tr th - { - border: 1px solid gray; - } - - table.simplelist tr td - { - border: none !important; - } - - /* Misc */ - span.highlight - { - font-weight: bold; - } - } - -/*============================================================================= - Images -=============================================================================*/ - - span.inlinemediaobject img - { - vertical-align: middle; - } - -/*============================================================================== - Super and Subscript: style so that line spacing isn't effected, see - http://www.adobe.com/cfusion/communityengine/index.cfm?event=showdetails&productId=1&postId=5341 -==============================================================================*/ - -sup, -sub { - height: 0; - line-height: 1; - vertical-align: baseline; - position: relative; - -} - -/* For internet explorer: */ - -* html sup, -* html sub { - vertical-align: bottom; -} - -sup { - bottom: 1ex; -} - -sub { - top: .5ex; -} - diff --git a/doc/images/alert.png b/doc/images/alert.png deleted file mode 100644 index b4645bc7..00000000 Binary files a/doc/images/alert.png and /dev/null differ diff --git a/doc/images/blank.png b/doc/images/blank.png deleted file mode 100644 index 764bf4f0..00000000 Binary files a/doc/images/blank.png and /dev/null differ diff --git a/doc/images/callouts/1.png b/doc/images/callouts/1.png deleted file mode 100644 index 6003ad3a..00000000 Binary files a/doc/images/callouts/1.png and /dev/null differ diff --git a/doc/images/callouts/1.svg b/doc/images/callouts/1.svg deleted file mode 100644 index e2e87dc5..00000000 --- a/doc/images/callouts/1.svg +++ /dev/null @@ -1,15 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/10.png b/doc/images/callouts/10.png deleted file mode 100644 index 0426f516..00000000 Binary files a/doc/images/callouts/10.png and /dev/null differ diff --git a/doc/images/callouts/10.svg b/doc/images/callouts/10.svg deleted file mode 100644 index 4740f587..00000000 --- a/doc/images/callouts/10.svg +++ /dev/null @@ -1,18 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/11.png b/doc/images/callouts/11.png deleted file mode 100644 index 821afc4f..00000000 Binary files a/doc/images/callouts/11.png and /dev/null differ diff --git a/doc/images/callouts/11.svg b/doc/images/callouts/11.svg deleted file mode 100644 index 09a0b2cf..00000000 --- a/doc/images/callouts/11.svg +++ /dev/null @@ -1,16 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/12.png b/doc/images/callouts/12.png deleted file mode 100644 index 7cec7272..00000000 Binary files a/doc/images/callouts/12.png and /dev/null differ diff --git a/doc/images/callouts/12.svg b/doc/images/callouts/12.svg deleted file mode 100644 index 9794044c..00000000 --- a/doc/images/callouts/12.svg +++ /dev/null @@ -1,18 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/13.png b/doc/images/callouts/13.png deleted file mode 100644 index 5b41e02a..00000000 Binary files a/doc/images/callouts/13.png and /dev/null differ diff --git a/doc/images/callouts/13.svg b/doc/images/callouts/13.svg deleted file mode 100644 index 64268bb4..00000000 --- a/doc/images/callouts/13.svg +++ /dev/null @@ -1,20 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/14.png b/doc/images/callouts/14.png deleted file mode 100644 index de5bdbd3..00000000 Binary files a/doc/images/callouts/14.png and /dev/null differ diff --git a/doc/images/callouts/14.svg b/doc/images/callouts/14.svg deleted file mode 100644 index 469aa974..00000000 --- a/doc/images/callouts/14.svg +++ /dev/null @@ -1,17 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/15.png b/doc/images/callouts/15.png deleted file mode 100644 index 3fd6ac38..00000000 Binary files a/doc/images/callouts/15.png and /dev/null differ diff --git a/doc/images/callouts/15.svg b/doc/images/callouts/15.svg deleted file mode 100644 index 8202233e..00000000 --- a/doc/images/callouts/15.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/16.svg b/doc/images/callouts/16.svg deleted file mode 100644 index 01d6bf81..00000000 --- a/doc/images/callouts/16.svg +++ /dev/null @@ -1,20 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/17.svg b/doc/images/callouts/17.svg deleted file mode 100644 index 0a04c556..00000000 --- a/doc/images/callouts/17.svg +++ /dev/null @@ -1,17 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/18.svg b/doc/images/callouts/18.svg deleted file mode 100644 index 1cb891b3..00000000 --- a/doc/images/callouts/18.svg +++ /dev/null @@ -1,21 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/19.svg b/doc/images/callouts/19.svg deleted file mode 100644 index e6fbb179..00000000 --- a/doc/images/callouts/19.svg +++ /dev/null @@ -1,20 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/2.png b/doc/images/callouts/2.png deleted file mode 100644 index f7c15788..00000000 Binary files a/doc/images/callouts/2.png and /dev/null differ diff --git a/doc/images/callouts/2.svg b/doc/images/callouts/2.svg deleted file mode 100644 index 07d03395..00000000 --- a/doc/images/callouts/2.svg +++ /dev/null @@ -1,17 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/20.svg b/doc/images/callouts/20.svg deleted file mode 100644 index ccbfd403..00000000 --- a/doc/images/callouts/20.svg +++ /dev/null @@ -1,20 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/21.svg b/doc/images/callouts/21.svg deleted file mode 100644 index 93ec53fd..00000000 --- a/doc/images/callouts/21.svg +++ /dev/null @@ -1,18 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/22.svg b/doc/images/callouts/22.svg deleted file mode 100644 index f48c5f3f..00000000 --- a/doc/images/callouts/22.svg +++ /dev/null @@ -1,20 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/23.svg b/doc/images/callouts/23.svg deleted file mode 100644 index 66242129..00000000 --- a/doc/images/callouts/23.svg +++ /dev/null @@ -1,22 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/24.svg b/doc/images/callouts/24.svg deleted file mode 100644 index a3d55253..00000000 --- a/doc/images/callouts/24.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/25.svg b/doc/images/callouts/25.svg deleted file mode 100644 index 56614a97..00000000 --- a/doc/images/callouts/25.svg +++ /dev/null @@ -1,21 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/26.svg b/doc/images/callouts/26.svg deleted file mode 100644 index 56faeaca..00000000 --- a/doc/images/callouts/26.svg +++ /dev/null @@ -1,22 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/27.svg b/doc/images/callouts/27.svg deleted file mode 100644 index a75c8121..00000000 --- a/doc/images/callouts/27.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/28.svg b/doc/images/callouts/28.svg deleted file mode 100644 index 7f8cf1a3..00000000 --- a/doc/images/callouts/28.svg +++ /dev/null @@ -1,23 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/29.svg b/doc/images/callouts/29.svg deleted file mode 100644 index cb63adf1..00000000 --- a/doc/images/callouts/29.svg +++ /dev/null @@ -1,22 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/3.png b/doc/images/callouts/3.png deleted file mode 100644 index 3ff0a939..00000000 Binary files a/doc/images/callouts/3.png and /dev/null differ diff --git a/doc/images/callouts/3.svg b/doc/images/callouts/3.svg deleted file mode 100644 index 918be806..00000000 --- a/doc/images/callouts/3.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/30.svg b/doc/images/callouts/30.svg deleted file mode 100644 index dc43ba1e..00000000 --- a/doc/images/callouts/30.svg +++ /dev/null @@ -1,22 +0,0 @@ - - - - -]> - - - - - - - - - diff --git a/doc/images/callouts/4.png b/doc/images/callouts/4.png deleted file mode 100644 index 6aa29fc0..00000000 Binary files a/doc/images/callouts/4.png and /dev/null differ diff --git a/doc/images/callouts/4.svg b/doc/images/callouts/4.svg deleted file mode 100644 index 8eb6a53b..00000000 --- a/doc/images/callouts/4.svg +++ /dev/null @@ -1,16 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/5.png b/doc/images/callouts/5.png deleted file mode 100644 index 36e78586..00000000 Binary files a/doc/images/callouts/5.png and /dev/null differ diff --git a/doc/images/callouts/5.svg b/doc/images/callouts/5.svg deleted file mode 100644 index ca7a9f22..00000000 --- a/doc/images/callouts/5.svg +++ /dev/null @@ -1,18 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/6.png b/doc/images/callouts/6.png deleted file mode 100644 index c943676b..00000000 Binary files a/doc/images/callouts/6.png and /dev/null differ diff --git a/doc/images/callouts/6.svg b/doc/images/callouts/6.svg deleted file mode 100644 index 783a0b9d..00000000 --- a/doc/images/callouts/6.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/7.png b/doc/images/callouts/7.png deleted file mode 100644 index 20940de3..00000000 Binary files a/doc/images/callouts/7.png and /dev/null differ diff --git a/doc/images/callouts/7.svg b/doc/images/callouts/7.svg deleted file mode 100644 index 59b3714b..00000000 --- a/doc/images/callouts/7.svg +++ /dev/null @@ -1,16 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/8.png b/doc/images/callouts/8.png deleted file mode 100644 index d8e34d4a..00000000 Binary files a/doc/images/callouts/8.png and /dev/null differ diff --git a/doc/images/callouts/8.svg b/doc/images/callouts/8.svg deleted file mode 100644 index c1803a3c..00000000 --- a/doc/images/callouts/8.svg +++ /dev/null @@ -1,20 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/callouts/9.png b/doc/images/callouts/9.png deleted file mode 100644 index abe63607..00000000 Binary files a/doc/images/callouts/9.png and /dev/null differ diff --git a/doc/images/callouts/9.svg b/doc/images/callouts/9.svg deleted file mode 100644 index bc149d3c..00000000 --- a/doc/images/callouts/9.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - -]> - - - - - - - - diff --git a/doc/images/caution.png b/doc/images/caution.png deleted file mode 100644 index 5b7809ca..00000000 Binary files a/doc/images/caution.png and /dev/null differ diff --git a/doc/images/caution.svg b/doc/images/caution.svg deleted file mode 100644 index 4bd586a0..00000000 --- a/doc/images/caution.svg +++ /dev/null @@ -1,68 +0,0 @@ - - - - - - Attenzione - - - - pulsante - - - - - Open Clip Art Library - - - - - Architetto Francesco Rollandin - - - - - Architetto Francesco Rollandin - - - - image/svg+xml - - - en - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/doc/images/draft.png b/doc/images/draft.png deleted file mode 100644 index 0084708c..00000000 Binary files a/doc/images/draft.png and /dev/null differ diff --git a/doc/images/home.png b/doc/images/home.png deleted file mode 100644 index 5584aacb..00000000 Binary files a/doc/images/home.png and /dev/null differ diff --git a/doc/images/home.svg b/doc/images/home.svg deleted file mode 100644 index e803a317..00000000 --- a/doc/images/home.svg +++ /dev/null @@ -1,26 +0,0 @@ - - - - - - - - -]> - - - - - - - - - - - - - - diff --git a/doc/images/important.png b/doc/images/important.png deleted file mode 100644 index 12c90f60..00000000 Binary files a/doc/images/important.png and /dev/null differ diff --git a/doc/images/important.svg b/doc/images/important.svg deleted file mode 100644 index dd84f3fe..00000000 --- a/doc/images/important.svg +++ /dev/null @@ -1,25 +0,0 @@ - - - - - - - - -]> - - - - - - - - - - - - - - - diff --git a/doc/images/next.png b/doc/images/next.png deleted file mode 100644 index 59800b4e..00000000 Binary files a/doc/images/next.png and /dev/null differ diff --git a/doc/images/next.svg b/doc/images/next.svg deleted file mode 100644 index 75fa83ed..00000000 --- a/doc/images/next.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - - - -]> - - - - - - - - - - - diff --git a/doc/images/next_disabled.png b/doc/images/next_disabled.png deleted file mode 100644 index 10a8c59d..00000000 Binary files a/doc/images/next_disabled.png and /dev/null differ diff --git a/doc/images/note.png b/doc/images/note.png deleted file mode 100644 index d0c3c645..00000000 Binary files a/doc/images/note.png and /dev/null differ diff --git a/doc/images/note.svg b/doc/images/note.svg deleted file mode 100644 index 648299d2..00000000 --- a/doc/images/note.svg +++ /dev/null @@ -1,33 +0,0 @@ - - - - - - - - - - - - -]> - - - - - - - - - - - - - - - - - - - diff --git a/doc/images/prev.png b/doc/images/prev.png deleted file mode 100644 index d88a40f9..00000000 Binary files a/doc/images/prev.png and /dev/null differ diff --git a/doc/images/prev.svg b/doc/images/prev.svg deleted file mode 100644 index 6d88ffdd..00000000 --- a/doc/images/prev.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - - - -]> - - - - - - - - - - - diff --git a/doc/images/prev_disabled.png b/doc/images/prev_disabled.png deleted file mode 100644 index ab3c17e0..00000000 Binary files a/doc/images/prev_disabled.png and /dev/null differ diff --git a/doc/images/smiley.png b/doc/images/smiley.png deleted file mode 100644 index 30a77f71..00000000 Binary files a/doc/images/smiley.png and /dev/null differ diff --git a/doc/images/tip.png b/doc/images/tip.png deleted file mode 100644 index 5c4aab3b..00000000 Binary files a/doc/images/tip.png and /dev/null differ diff --git a/doc/images/tip.svg b/doc/images/tip.svg deleted file mode 100644 index cd437a5e..00000000 --- a/doc/images/tip.svg +++ /dev/null @@ -1,84 +0,0 @@ - - - - - - lamp - - - - office - - lamp - - - - - Open Clip Art Library - - - - - Sergio Luiz Araujo Silva - - - - - Public Domain - - - set 2005 - image/svg+xml - - - en - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/doc/images/toc-blank.png b/doc/images/toc-blank.png deleted file mode 100644 index 6ffad17a..00000000 Binary files a/doc/images/toc-blank.png and /dev/null differ diff --git a/doc/images/toc-minus.png b/doc/images/toc-minus.png deleted file mode 100644 index abbb020c..00000000 Binary files a/doc/images/toc-minus.png and /dev/null differ diff --git a/doc/images/toc-plus.png b/doc/images/toc-plus.png deleted file mode 100644 index 941312ce..00000000 Binary files a/doc/images/toc-plus.png and /dev/null differ diff --git a/doc/images/up.png b/doc/images/up.png deleted file mode 100644 index 17d9c3ec..00000000 Binary files a/doc/images/up.png and /dev/null differ diff --git a/doc/images/up.svg b/doc/images/up.svg deleted file mode 100644 index d31aa9c8..00000000 --- a/doc/images/up.svg +++ /dev/null @@ -1,19 +0,0 @@ - - - - - - -]> - - - - - - - - - - - diff --git a/doc/images/up_disabled.png b/doc/images/up_disabled.png deleted file mode 100644 index e22bc871..00000000 Binary files a/doc/images/up_disabled.png and /dev/null differ diff --git a/doc/images/warning.png b/doc/images/warning.png deleted file mode 100644 index 1c33db8f..00000000 Binary files a/doc/images/warning.png and /dev/null differ diff --git a/doc/images/warning.svg b/doc/images/warning.svg deleted file mode 100644 index fc8d7484..00000000 --- a/doc/images/warning.svg +++ /dev/null @@ -1,23 +0,0 @@ - - - - - - - - -]> - - - - - - - - - - - - - diff --git a/doc/index.html b/doc/index.html deleted file mode 100644 index f3eac1d4..00000000 --- a/doc/index.html +++ /dev/null @@ -1,103 +0,0 @@ - - - -Chapter 1. Boost.Numeric.Odeint - - - - - - -
Next
-
-
-

-Chapter 1. Boost.Numeric.Odeint

-

-Karsten Ahnert -

-

-Mario Mulansky -

-

Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

-
-

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-

Table of Contents

-
-
Getting started
-
-
Overview
-
Usage, - Compilation, Headers
-
Short - Example
-
-
Tutorial
-
-
Harmonic - oscillator
-
Solar system
-
Chaotic - systems and Lyapunov exponents
-
Stiff systems
-
Special - topics
-
Using - Cuda and Thrust
-
Using - OpenCL and VexCL
-
All examples
-
References
-
-
odeint in detail
-
-
Steppers
-
Generation - functions
-
Integrate - functions
-
Iterators
-
State - types, algebras and operations
-
Using - boost::ref
-
Using - boost::range
-
-
Concepts
-
-
System
-
Symplectic - System
-
Simple - Symplectic System
-
Implicit - System
-
Stepper
-
Error Stepper
-
Controlled - Stepper
-
Dense - Output Stepper
-
State - Algebra Operations
-
State Wrapper
-
-
-
-
- - - -

Last revised: September 07, 2012 at 19:18:24 GMT

-
-
Next
- - diff --git a/doc/odeint/reference.html b/doc/odeint/reference.html deleted file mode 100644 index e773413a..00000000 --- a/doc/odeint/reference.html +++ /dev/null @@ -1,33 +0,0 @@ - - - -Reference - - - - - - - -
-
-
-PrevUpHome -
-

-Reference

- - - -
Copyright © 2009-2011 Karsten Ahnert - and Mario Mulansky

- Distributed under the Boost Software License, Version 1.0. (See accompanying - file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) -

-
-
-
-PrevUpHome -
- - diff --git a/doc/odeint_HTML.manifest b/doc/odeint_HTML.manifest deleted file mode 100644 index 007493e6..00000000 --- a/doc/odeint_HTML.manifest +++ /dev/null @@ -1,49 +0,0 @@ -index.html -boost_sandbox_numeric_odeint/getting_started.html -boost_sandbox_numeric_odeint/getting_started/overview.html -boost_sandbox_numeric_odeint/getting_started/usage__compilation__headers.html -boost_sandbox_numeric_odeint/getting_started/short_example.html -boost_sandbox_numeric_odeint/tutorial.html -boost_sandbox_numeric_odeint/tutorial/harmonic_oscillator.html -boost_sandbox_numeric_odeint/tutorial/solar_system.html -boost_sandbox_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html -boost_sandbox_numeric_odeint/tutorial/stiff_systems.html -boost_sandbox_numeric_odeint/tutorial/special_topics.html -boost_sandbox_numeric_odeint/tutorial/using_cuda_and_thrust.html -boost_sandbox_numeric_odeint/tutorial/all_examples.html -boost_sandbox_numeric_odeint/tutorial/references.html -boost_sandbox_numeric_odeint/odeint_in_detail.html -boost_sandbox_numeric_odeint/odeint_in_detail/steppers.html -boost_sandbox_numeric_odeint/odeint_in_detail/integrate_functions.html -boost_sandbox_numeric_odeint/odeint_in_detail/algebras_and_operations.html -boost_sandbox_numeric_odeint/odeint_in_detail/using_boost__ref.html -boost_sandbox_numeric_odeint/odeint_in_detail/using_boost__range.html -boost_sandbox_numeric_odeint/extend_odeint.html -boost_sandbox_numeric_odeint/extend_odeint/write_own_steppers.html -boost_sandbox_numeric_odeint/extend_odeint/adapt_your_own_containers.html -boost_sandbox_numeric_odeint/extend_odeint/adapt_your_own_operations.html -boost_sandbox_numeric_odeint/concepts.html -boost_sandbox_numeric_odeint/concepts/system.html -boost_sandbox_numeric_odeint/concepts/symplectic_system.html -boost_sandbox_numeric_odeint/concepts/simple_symplectic_system.html -boost_sandbox_numeric_odeint/concepts/implicit_system.html -boost_sandbox_numeric_odeint/concepts/observer.html -boost_sandbox_numeric_odeint/concepts/stepper.html -boost_sandbox_numeric_odeint/concepts/error_stepper.html -boost_sandbox_numeric_odeint/concepts/controlled_stepper.html -boost_sandbox_numeric_odeint/concepts/dense_output_stepper.html -boost_sandbox_numeric_odeint/concepts/state_algebra_operations.html -odeint/reference.html -boost_sandbox_numeric_odeint/old_concepts.html -boost_sandbox_numeric_odeint/old_concepts/basic_stepper.html -boost_sandbox_numeric_odeint/old_concepts/error_stepper.html -boost_sandbox_numeric_odeint/old_concepts/controlled_stepper.html -boost_sandbox_numeric_odeint/old_concepts/dense_ouput_stepper.html -boost_sandbox_numeric_odeint/old_concepts/size_adjusting_stepper.html -boost_sandbox_numeric_odeint/old_concepts/compositestepper.html -boost_sandbox_numeric_odeint/old_reference.html -boost_sandbox_numeric_odeint/old_reference/stepper_classes.html -boost_sandbox_numeric_odeint/old_reference/integration_functions.html -boost_sandbox_numeric_odeint/old_reference/algebras.html -boost_sandbox_numeric_odeint/old_reference/operations.html -boost_sandbox_numeric_odeint/old_reference/resizing.html diff --git a/doc/phase_lattice_2d_0000.jpg b/doc/phase_lattice_2d_0000.jpg deleted file mode 100644 index f98703eb..00000000 Binary files a/doc/phase_lattice_2d_0000.jpg and /dev/null differ diff --git a/doc/phase_lattice_2d_0100.jpg b/doc/phase_lattice_2d_0100.jpg deleted file mode 100644 index 78628393..00000000 Binary files a/doc/phase_lattice_2d_0100.jpg and /dev/null differ diff --git a/doc/phase_lattice_2d_1000.jpg b/doc/phase_lattice_2d_1000.jpg deleted file mode 100644 index e0ee037c..00000000 Binary files a/doc/phase_lattice_2d_1000.jpg and /dev/null differ diff --git a/doc/solar_system.jpg b/doc/solar_system.jpg deleted file mode 100644 index b89c3b3d..00000000 Binary files a/doc/solar_system.jpg and /dev/null differ diff --git a/libs/numeric/odeint/test/prepare_stepper_testing.hpp b/libs/numeric/odeint/test/prepare_stepper_testing.hpp index 9669df4e..fda010a5 100644 --- a/libs/numeric/odeint/test/prepare_stepper_testing.hpp +++ b/libs/numeric/odeint/test/prepare_stepper_testing.hpp @@ -17,7 +17,7 @@ #ifndef PREPARE_STEPPER_TESTING_HPP_ #define PREPARE_STEPPER_TESTING_HPP_ -#include +#include #include #include @@ -56,8 +56,8 @@ struct constant_system_fusion typedef mpl::vector < - mpl::vector< float , std::tr1::array< float , 1 > , std::tr1::array< float , 1 > , constant_system_standard , stepper_type > , - mpl::vector< float , std::tr1::array< float , 1 > , std::vector< float > , constant_system_standard , stepper_type > + mpl::vector< float , boost::array< float , 1 > , boost::array< float , 1 > , constant_system_standard , stepper_type > , + mpl::vector< float , boost::array< float , 1 > , std::vector< float > , constant_system_standard , stepper_type > > types_and_systems_matrix;