pBook‎ > ‎

### Solver:: Integrator:: C++

#### Bug in odeint in Boost below 1.58

posted Dec 2, 2014, 8:32 AM by Javad Taghia   [ updated Dec 2, 2014, 8:33 AM ]

## https://svn.boost.org/trac/boost/ticket/10487#comment:3Multiple definitions of boost::numeric::odeint::detail::* by multiple #include

Reported by: Owned by: guillaume.jacquenot@… karsten Boost 1.57.0 odeint Boost 1.56.0 Problem guillaume.jacquenot@…

### Description

Dear all,

I encountered a linker problem within the odeint code of library boost::numeric with boost 1.56 on Debian 32 bit with gcc 4.7. There (and I guess on many other compilers) will be multiple definitions of functions. The error is known not to occur with version 1.55.0.

When using several inclusions of the following files in various code file

#include <boost/numeric/odeint/stepper/euler.hpp> #include <boost/numeric/odeint/stepper/runge_kutta4.hpp> #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>

the linker produces an error of multiple definitions for boost::numeric::odeint::detail.

This linker problem will not occur if only one .cpp file #includes this file.

Possible explanation taken from bug number 7678:

Seemingly compiling each compile unit like foo.cpp, bar.cpp, etc. produces a definition into its individual .obj-file, namely foo.obj, bar.obj, etc.. Later the linker can not resolve which of the multiple definitions of "belongs", the one in foo.obj or the one in bar.obj, to take when binding the executable.

CMake log:

Linking CXX executable run_all_tests XXXXXXXXXXX/libXXXXXXXXXXX_static.a(FILE_YYYYYYYYYYYYY.cpp.o):(.bss+0x0): multiple definition of boost::numeric::odeint::detail::_2' src/FILE_ZZZZZZZZZZZ.cpp.o:(.bss+0x10): first defined here XXXXXXXXXXX/libXXXXXXXXXXX_static.a(FILE_YYYYYYYYYYYYY.cpp.o):(.bss+0x1): multiple definition of boost::numeric::odeint::detail::_1' src/FILE_ZZZZZZZZZZZ.cpp.o:(.bss+0x11): first defined here collect2: error: ld returned 1 exit status

The bug is similar to the one described here

https://svn.boost.org/trac/boost/ticket/7678

Guillaume Jacquenot

## Change History

### comment:1Changed 3 months ago by guillaume.jacquenot@…

• Summary changed from Multiple definitions of bool boost::numeric::odeint by multiple #include toMultiple definitions of boost::numeric::odeint::detail::* by multiple #include

### comment:2Changed 3 weeks ago by karsten

• Status changed from new to closed
• Resolution set to fixed

The fix will apear in Boost 1.58.

### comment:3Changed 1 second ago by telcom

The same issue exists in Visual Studio 2012 C++11, Boost 1.56. Possible temporary solution is, getting odeint-v2 from gitub and replace the following folders and file in the include path of numeric: ~/numeric/odeint :: folder odeint.hpp :: file

Tested and verified.

#### Chaotic Systems and odeint

posted Dec 2, 2014, 6:13 AM by Javad Taghia

### Chaotic systems and Lyapunov exponents

In this example we present application of odeint to investigation of 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 exist, 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 the 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 convergence to a fixed point is expected. 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 of 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 of lorenz() and lorenz_with_lyap() having different state types may be solved putting 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:

// explicitly choose range_algebra to override default choice of array_algebra
runge_kutta4< state_type , double , state_type , double , range_algebra > rk4;

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

Note Note that when using Boost.Range, we have to explicitly configure the stepper to use the range_algebra as otherwise odeint would automatically chose the array_algebra, which is incompatible with the usage of Boost.Range, because the original state_type is an array.

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

1. Initialize the perturbations. They are stored linearly behind the state of the Lorenz system. The perturbations are initialized such that ​ij = δ ​ij, where ​ij is the j-component of the i.-th perturbation and δ ​ij is the Kronecker symbol.
2. Integrate 100 steps of the full system with perturbations
3. Orthonormalize the perturbation using Gram-Schmidt orthonormalization algorithm.
4. Repeat step 2 and 3. Every 10000 steps write the current Lyapunov exponent.
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;
}
}

#### Example Kin Sim

posted Dec 1, 2014, 10:08 PM by Javad Taghia   [ updated Dec 2, 2014, 8:26 PM ]

 //============================================================================//============================================================================#include #include #include #include #include #include typedef std::vector state_type;const double len = 1.7;// the rhs of x' = f(x)void harmonic_oscillator(const state_type &x, state_type &dxdt, const double t) { dxdt[0] =  x[3] * cos(x[2]); dxdt[1] =  x[3] * sin(x[2]); dxdt[2] = x[3] / len * tan(x[4]); dxdt[3] = 0; dxdt[4] = 0; }template std::vector createRange(vectype start, vectype end, vectype step ){    std::vector result;    for (; start < end; start += step) {        result.push_back(start);    }    return result;}int main() { state_type x(5); x[0] = 0.0; // xt x[1] = 0.0; // yt x[2] = 0.0; // tHeta x[3] = 1; // vf x[4] = 0; // dElta_t std::vector t = createRange(0, 10, 0.01); //for (double i=0; i<=10; i+=0.1) //t.push_back(i); std::cout << x[0] << "  " << x[1] << "  " << x[2] << std::endl; double t_old = 0; for (auto i: t) { size_t steps = boost::numeric::odeint:: integrate(harmonic_oscillator, x, t_old, i, 0.001); t_old = i; std::cout << x[0] << "  " << x[1] << "  " << x[2] << std::endl; } return 0;}with array#include #include #include typedef  boost::array state_type;const double len = 1.7;int main() { ashacontrol::stateModel model(len); state_type x; x[0] = 0.0; // xt x[1] = 0.0; // yt x[2] = 0.0; // tHeta x[3] = 1; // vf x[4] = 0; // dElta_t std::vector t = ashacontrol::createRange(0, 10, 0.01); std::cout << x[0] << "  " << x[1] << "  " << x[2] << std::endl; double t_old = 0; for (auto i: t) { size_t steps = boost::numeric::odeint:: integrate(model, x, t_old, i, 0.001); t_old = i; std::cout << x[0] << "  " << x[1] << "  " << x[2] << std::endl; }   getchar(); ashacontrol::CStopwatch timer1; splab::Timing splabtimer2; timer1.start(); splabtimer2.start(); std::string line; boost::regex pat( "^Subject: (Re: |Aw: )*(.*)" ); while (std::cin) { std::cout<<"hi"<

#### odeint2

posted Nov 29, 2014, 6:32 AM by Javad Taghia

odeint is a header-only library, no linking against pre-compiled code is required. It can be included 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;


It is also possible to include only parts of the library. This is the recommended way since it saves a lot of compilation time.

• #include <boost/numeric/odeint/stepper/XYZ.hpp> - the include path for all steppers, XYZ is a placeholder for a stepper.
• #include <boost/numeric/odeint/algebra/XYZ.hpp> - all algebras.
• #include <boost/numeric/odeint/util/XYZ.hpp> - the utility functions like is_resizeable,same_size, or resize.
• #include <boost/numeric/odeint/integrate/XYZ.hpp> - the integrate routines.
• #include <boost/numeric/odeint/iterator/XYZ.hpp> - the range and iterator functions.
• #include <boost/numeric/odeint/external/XYZ.hpp> - any binders to external libraries.

### 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 exampleboost::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_kutta54_cash_karp 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 also possible to represent the ode system as a class. The rhs must then be implemented as a functor - a class with an overloaded function call 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 );


In order to observe the solution during the integration steps 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. 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: harmonic_oscillator.cpp

#### odeint

posted Nov 29, 2014, 6:29 AM by Javad Taghia

### Overview

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

Stepper

System

configurable

No

No

Yes

Multistep method

adams_moulton

Stepper

System

configurable

No

No

Yes

Multistep method

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 eulerand 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 withBoost.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.uBLASvectors 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.uBLASvectors as state types.

Symplectic Euler

symplectic_euler

Stepper

1

No

No

No

Basic symplectic solver for separable Hamiltonian system

Symplectic RKN McLachlan

symplectic_rkn_sb3a_mclachlan

Stepper

4

No

No

No

Symplectic solver for separable Hamiltonian system with 6 stages and order 4.

Symplectic RKN McLachlan

symplectic_rkn_sb3a_m4_mclachlan

Stepper

4

No

No

No

Symplectic solver with 5 stages and order 4, can be used with arbitrary precision types.

Velocity Verlet

velocity_verlet

Stepper

Second Order System

1

No

No

Yes

Velocity verlet method suitable for molecular dynamics simulation.

#### Boost

posted Nov 29, 2014, 5:24 AM by Javad Taghia   [ updated Nov 29, 2014, 5:24 AM ]

 In C++11  a simple lambda function wrapping the call to your member methodClass c; auto f = [&c]( const state_type & x , state_type &dxdt , double t ) { c.system_func( x , dxdt , t ); }; integrate( f , x0 , t0 , t1 , dt );std::bind might also work, but then you have to take care if values are passed by reference of by value.In C++03 you need to write a simple wrapper around your class methodstruct wrapper { Class &c; wrapper( Class &c_ ) : c( c_ ) { } template< typename State , typename Time > void operator()( State const &x , State &dxdt , Time t ) const { c.system_func( x , dxdt , t ); } }; // ... integrate( wrapper( c ) , x0 , t0 , t1 , dt );(Boost.Bind will not work correctly with more than two arguments).#include #include #include using namespace std; using namespace boost::numeric::odeint; const double sigma = 10.0; const double R = 28.0; const double b = 8.0 / 3.0; typedef boost::array< double , 3 > state_type; void lorenz( const state_type &x , 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]; } void write_lorenz( const state_type &x , const double t ) { cout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl; } int main(int argc, char **argv) { state_type x = { 10.0 , 1.0 , 1.0 }; // initial conditions integrate( lorenz , x , 0.0 , 25.0 , 0.1 , write_lorenz ); }And regarding the system, you can provide anything where the following expression is valid:sys( x , dxdt , t ) // returning void`Check the user's guide (and more examples) online.

#### C++11

posted Nov 29, 2014, 5:22 AM by Javad Taghia   [ updated Nov 29, 2014, 5:22 AM ]

https://github.com/michaeloriordan/integrator

 #include "Integrator.hpp" #include "VecDef.hpp" #include #include void fillVec(vec_t& vec, std::function f, double dx) { const int N = vec.size(); for (int i = 0; i < N; ++i) { vec[i] = f(i * dx); } } int main() { double a = 0; double b = 2; int N = 1000; double dx = (b - a) / N; vec_t vec(N, 0.0); auto f = [](double x){return x*x;}; fillVec(vec, f, dx); double sum = integrator::integrate(vec, dx); std::cout << sum << std::endl; return 0; }

 #ifndef __Integrator_ #define __Integrator_ #include #include #include #include #include namespace si_lib { /** The calendar queue class stores a (time,event_name,recurrence_interval) tuple. It is essentially a priority queue whose top priority is always increasing. */ class CalendarQueue { private: /** This class stores an event. An event is represented by a floating-point time and a string identifying what the event is. The event may recur at a regular interval. */ class DiscreteEvent { public: double t; std::string event; double recur_int; DiscreteEvent (double t, const char event[], double recur_int) : t(t), event(event), recur_int(recur_int) {} DiscreteEvent (double t, const std::string &event, double recur_int) : t(t), event(event), recur_int(recur_int) {} bool operator< (const DiscreteEvent& a) const { return t>a.t; } }; typedef std::priority_queue > deq_type; deq_type deq; public: ///Schedule \a event at \a t void insert(double t, const char event[], double recur_int=0){ assert(recur_int>=0); //assert(empty() || t>=current_time()); deq.push(DiscreteEvent(t,event,recur_int)); } ///Schedule \a event at \a t void insert(double t, const std::string &event, double recur_int=0){ assert(recur_int>=0); //assert(empty() || t>=current_time()); deq.push(DiscreteEvent(t,event,recur_int)); } ///Time the top element of the calendar queue is scheduled for double current_time() const { return deq.top().t; } ///Event indicated by the top element of the calendar queue std::string current_event() const { return deq.top().event; } ///Remove the top element of the calendar queue, and cast it into oblivion void pop() { if(deq.empty()) return; if(deq.top().recur_int>0) insert(current_time()+deq.top().recur_int, deq.top().event, deq.top().recur_int); deq.pop(); } ///Take the top element of the queue and insert a copy of it \a dt ///into the future double reschedule_top(double dt) { deq.push(DiscreteEvent(deq.top().t+dt, deq.top().event, deq.top().recur_int)); return current_time()+dt; } ///Clear all events from the calendar queue void clear() { deq=deq_type(); }; ///Is the calendar queue empty? bool empty() const { return deq.empty(); } }; /** The ArrayState class wraps std::array providing algebra operations used by the Integrator. */ template class ArrayState : public std::array { public: ArrayState() {} ArrayState(const std::vector &init){ for(unsigned int i=0;i::size();++i) std::array::operator[](i)=init[i]; } ArrayState& operator+=(const ArrayState &other) { for(unsigned int i=0;i::size();++i) std::array::operator[](i)+=other[i]; return *this; } }; template ArrayState operator+( const ArrayState &a , double b ){ ArrayState result(a); for(T& i: result) i+=b; return result; } template ArrayState operator+( double a, const ArrayState &b ){ return b+a; } template ArrayState operator+( const ArrayState &a, const ArrayState &b ){ ArrayState result(a); for(unsigned int i=0;i ArrayState operator*( const ArrayState &a, double b ){ ArrayState result(a); for(T& i: result) i*=b; return result; } template ArrayState operator*( double a, const ArrayState &b ){ return b*a; } template ArrayState operator/( const ArrayState &a, double b){ ArrayState result(a); for(T& i: result) i/=b; return result; } template double abs( const ArrayState &a ){ double result=0; for(const T& i: a) result+=abs(i); return result; } template double abs( const ArrayState &a ){ double result=0; for(const double& i: a) result+=std::abs(i); return result; } template std::ostream& operator<<(std::ostream &out, const ArrayState &a){ for(typename ArrayState::const_iterator i=a.begin();i!=a.end();++i) out<<*i<<" "; return out; } /** The Integrator class supplies a simple Euler's method integrator with an adaptive step-size scheme */ template class Integrator { protected: T stateval; typedef std::function dx_type; dx_type dx; double dtval, dtmin, dtmax; double t; int goodsteps; int stepcount; public: Integrator(const T &stateval, dx_type dx, double dtmin, double dtmax); double time() const; T& state(); double dt() const; void dt(double h); void step(); int steps() const; }; template Integrator::Integrator(const T &stateval, dx_type dx, double dtmin, double dtmax) : stateval(stateval), dx(dx), dtmin(dtmin), dtmax(dtmax) { assert(dtmin>0); assert(dtmax>0); assert(dtmin<=dtmax); dtval=dtmin; stepcount=0; t=0; } template double Integrator::time() const { return t; } template T& Integrator::state() { return stateval; } template double Integrator::dt() const { return dtval; } template void Integrator::dt(double h) { assert(h>0); assert(h>=dtmin); assert(h<=dtmax); dtval=h; } template int Integrator::steps() const { return stepcount; } template void Integrator::step() { T e1, e2; ++stepcount; dx(stateval , e1, t); dx(stateval+e1*dtval/2, e2, t+dtval/2.); double abs_e1=abs(stateval+e1*dtval); double abs_e2=abs(stateval+e1*dtval/2+e2*dtval/2); if( (std::abs(abs_e1-abs_e2)/(std::abs(abs_e1+abs_e2)/2))>0.05 ){ stateval+=e1*dtval/2.; t+=dtval/2.; dtval/=2.; if(dtval=8) dtval*=2; if(dtval>dtmax) dtval=dtmax; } /** The EventIntegrator class supplies a simple Euler's method integrator with an adaptive step-size scheme. It also allows for discrete events. It approaches these events with an exponentially-decreasing stepsize and moves away from them in the same fashion. */ template class EventIntegrator : public Integrator { private: CalendarQueue calq; bool at_event; std::string at_event_name; public: EventIntegrator(const T &stateval, typename Integrator::dx_type dx, double dtmin, double dtmax) : Integrator(stateval, dx, dtmin, dtmax) { at_event=false; at_event_name=""; } void insert_event(double t, const std::string &event, double recur_int); void insert_event(double t, const char event[], double recur_int); bool is_event() const; std::string event() const; void step(); }; template void EventIntegrator::insert_event(double t, const std::string &event, double recur_int=0){ calq.insert(t,event,recur_int); } template void EventIntegrator::insert_event(double t, const char event[], double recur_int=0){ calq.insert(t,event,recur_int); } template bool EventIntegrator::is_event() const { return at_event; } template std::string EventIntegrator::event() const { if( is_event() ) return calq.current_event(); else return ""; } template void EventIntegrator::step() { T e1, e2; if(at_event && calq.current_time()==Integrator::t){ at_event_name=calq.current_event(); calq.pop(); return; } at_event=false; ++Integrator::stepcount; if(!calq.empty()){ while( Integrator::dtval > Integrator::dtmin && Integrator::t + Integrator::dtval > calq.current_time() ) { Integrator::dtval/=2.; } if(Integrator::dtval < Integrator::dtmin) Integrator::dtval = Integrator::dtmin; if(Integrator::dtval == Integrator::dtmin && Integrator::t + Integrator::dtval > calq.current_time()){ Integrator::dtval=calq.current_time()-Integrator::t; Integrator::dx(Integrator::stateval, e1, Integrator::t); Integrator::stateval+=e1*Integrator::dtval; Integrator::t=calq.current_time(); Integrator::dtval = Integrator::dtmin; at_event=true; at_event_name=calq.current_event(); calq.pop(); return; } } Integrator::dx(Integrator::stateval , e1, Integrator::t); Integrator::dx(Integrator::stateval + e1*Integrator::dtval/2, e2, Integrator::t+Integrator::dtval/2.); double abs_e1=abs(Integrator::stateval + e1*Integrator::dtval); double abs_e2=abs(Integrator::stateval + e1*Integrator::dtval/2 + e2*Integrator::dtval/2); if( (std::abs(abs_e1-abs_e2)/(std::abs(abs_e1+abs_e2)/2))>0.05 ){ Integrator::stateval+=e1*Integrator::dtval/2.; Integrator::t+=Integrator::dtval/2.; Integrator::dtval/=2.; if(Integrator::dtval::dtmin) Integrator::dtval=Integrator::dtmin; Integrator::goodsteps=0; } else { Integrator::stateval+=e1*Integrator::dtval; Integrator::t+=Integrator::dtval; ++Integrator::goodsteps; } if(Integrator::dtval::dtmax && Integrator::goodsteps>=8) Integrator::dtval*=2; if(Integrator::dtval>Integrator::dtmax) Integrator::dtval=Integrator::dtmax; } } #endif

1-7 of 7