State type boost:fusion in boost:odeint

Question

I want to use boost::odeint to solve differential equations for different collections of variables- say of std::vector type- in parallel. One solution would of course be to combine all variables into a large vector which is then employed as state variable.

However, I would prefer a more elegant solution like employing boost::fusion as a state type which then holds the different vectors. As far as I understood from a posting and the implementation for a related problem, there is -in principle- no obstacle in doing so. I only miss a few hints for the concrete implementation- in particular concerning the right specification of

algebra, operations and resizing

required for the creation for example of an error stepper. Which of the existing implementations- e.g. odeint::fusion_algebra- can be used directly and what remains to done in this case?

Answer 1

Using boost::fusion compile time containers if fine with the fusion_algebra and the default_operations, as long as each of its elements support

• Multiplication with a scalar (*)
• Addition and subtraction (+,-)
• Resizing

This is the case for the elementary floating point types, like double, float, or even std::complex. This is also the case for types supporting expression templates like the vector and matrix types from the MTL, boost::ublas, or vexcl and viennacl. But it is not possible for std::vector, since vectors do not implement or overload the appropriate +*-/ operators. In this case you have three possibilities

1. Choose an vector type supporting expression templates
2. Implement a nested algebra which iterates the fusion sequence with fusion's for_each and calls the correct for_each from the sub algebra
3. Implement a custom operation which iterates over the elements of the fusion sequence.

Note 1: Resizing usually means, that if your type needs resizing you specialize is_resizable to a compile time true, for example via

template<>
is_resizable< YourType > : boost::true_type { };

and specialize resize_impl<> and same_size<>. Have a look at the default implementations. This is really ease and should be only a one-liner.

Note 2: If you also need step size control and you choose to take a vector type with expresson templates the operator / and the function max must exist and they must perform element-wise division and element-wise max. This might be tricky to implement. If you need some help feel free to contact us again.