Search code examples
c++boostodeint

State type boost:fusion in boost:odeint

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?

Solution

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