AdamsBashforth_fwd.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_INTEGRATORS_ADAMSBASHFORTH_FWD_HPP_
#define MEDUSA_BITS_INTEGRATORS_ADAMSBASHFORTH_FWD_HPP_
 
#include <medusa/Config.hpp>
#include <Eigen/Core>
#include <medusa/bits/utils/numutils.hpp>
#include "RKExplicit_fwd.hpp"
#include <medusa/bits/utils/assert.hpp>
 
namespace mm {
 
template <typename Scalar, int num_steps>
class AdamsBashforth {
  public:
    typedef Scalar scalar_t;  
    enum {
        steps = num_steps  
    };
  private:
    Eigen::Matrix<scalar_t, steps, 1> b_;  
 
  public:
    AdamsBashforth(const Eigen::Matrix<scalar_t, steps, 1>& b) : b_(b) {
        assert_msg(std::abs(b_.sum() - 1.0) < 1e-15, "AB method coefficients do not sum up to 1, "
                                                     "but %.16f instead.", b_.sum());
    }
 
    template <typename func_t, typename initial_method_t = RKExplicit<Scalar, num_steps>>
    class Integrator {
        const initial_method_t initial_method_;  
        const AdamsBashforth<scalar_t, num_steps> method_;  
        scalar_t t0_;  
        scalar_t dt_;  
        int max_steps;  
        Eigen::VectorXd y0_;  
        const func_t& func_;  
 
      public:
        Integrator(const AdamsBashforth<scalar_t, num_steps>& method,
                   const initial_method_t& initial_method,
                   const func_t& func, scalar_t t0, scalar_t t_max, scalar_t dt, Eigen::VectorXd y0)
                : initial_method_(initial_method), method_(method), t0_(t0),
                  dt_(dt), max_steps(iceil((t_max - t0) / dt)), y0_(std::move(y0)),
                  func_(func) {}
 
        class IterationStep : public std::iterator<
                std::forward_iterator_tag,   // iterator_category
                IterationStep,               // value_type
                int,                         // difference_type
                const IterationStep*,        // pointer
                IterationStep                // reference
        > {
            const Integrator& integrator;  
            scalar_t t;  
            Eigen::Matrix<scalar_t, Eigen::Dynamic, steps> last_ys;  
            int cur_step;  
 
          public:
            explicit IterationStep(const Integrator& integrator)
                    : integrator(integrator), t(integrator.t0_),
                      last_ys(integrator.y0_.size(), static_cast<int>(steps)), cur_step(0) {
                last_ys.col(steps-1) = integrator.y0_;
            }
 
            IterationStep(const Integrator& integrator, int) :
                    integrator(integrator), t(), last_ys(), cur_step(integrator.max_steps + 1) {}
 
            IterationStep& operator++();
 
            IterationStep operator++(int) {
                IterationStep retval = *this;
                ++(*this);
                return retval;
            }
 
            bool operator==(const IterationStep& other) const { return cur_step == other.cur_step; }
 
            bool operator!=(const IterationStep& other) const { return !(*this == other); }
 
            IterationStep& operator*() { return *this; }
 
            const IterationStep& operator*() const { return *this; }
 
            explicit operator bool() {
                return cur_step <= integrator.max_steps;
            }
 
            bool is_last() const {
                return cur_step == integrator.max_steps;
            }
 
            scalar_t time() const { return t; }
 
            typename Eigen::Matrix<scalar_t, Eigen::Dynamic, steps>::ConstColXpr value() const {
                return last_ys.col(steps-1);
            }
 
            scalar_t& time() { return t; }
 
            typename Eigen::Matrix<scalar_t, Eigen::Dynamic, steps>::ColXpr value() {
                return last_ys.col(steps-1);
            }
 
            friend std::ostream& operator<<(std::ostream& os, const IterationStep& step) {
                return os << "t = " << step.t << "; y = " << step.last_ys.col(steps-1);
            }
        };
 
        typedef IterationStep iterator;  
        typedef IterationStep const_iterator;  
 
        iterator begin() { return IterationStep(*this); }
 
        const_iterator cbegin() { return IterationStep(*this); }
 
        iterator end() { return IterationStep(*this, 0); }
 
        const_iterator cend() { return IterationStep(*this, 0); }
 
        friend std::ostream& operator<<(std::ostream& os, const Integrator& integrator) {
            return os << "Multistep integrator using " << integrator.method_ << " method from "
                      << "time " << integrator.t0_ << " with step " << integrator.dt_ << " for "
                      << integrator.max_steps << " steps.";
        }
    };
 
    template <typename func_t>
    Integrator<func_t> solve(const func_t& func, scalar_t t0, scalar_t t_max, scalar_t dt,
                             const Eigen::VectorXd& y0) const {
        return Integrator<func_t, RKExplicit<scalar_t, num_steps>>(
                *this, integrators::Explicit::of_order<num_steps>(), func, t0, t_max, dt, y0);
    }
 
    template <typename initial_method_t, typename func_t>
    Integrator<func_t, initial_method_t> solve(const initial_method_t& method, const func_t& func,
                                               scalar_t t0, scalar_t t_max, scalar_t dt,
                                               const Eigen::VectorXd& y0) const {
        return Integrator<func_t, initial_method_t>(*this, method, func, t0, t_max, dt, y0);
    }
 
  private:
    template <typename func_t>
    Eigen::VectorXd step(const func_t& func, scalar_t t,
                         const Eigen::Matrix<scalar_t, Eigen::Dynamic, steps>& last_ys,
                         scalar_t dt) const {
        Eigen::VectorXd result = last_ys.col(steps-1);
        for (int i = 0; i < steps; ++i) {
            result += dt * b_(i) * func(t - i*dt, last_ys.col(i));
        }
        return result;
    }
 
  public:
    template <typename scalar_t, int stages>
    friend std::ostream& operator<<(std::ostream& os, const AdamsBashforth<scalar_t, stages>&);
};
 
namespace integrators {
 
namespace ExplicitMultistep {
 
template <class scalar_t = double>
static AdamsBashforth<scalar_t, 1> AB1() {
    Eigen::Matrix<scalar_t, 1, 1> b(1); b << 1.0;
    return {b};
}
 
template <class scalar_t = double>
static AdamsBashforth<scalar_t, 2> AB2() {
    Eigen::Matrix<scalar_t, 2, 1> b(2); b << -1./2., 3./2.;
    return {b};
}
 
template <class scalar_t = double>
static AdamsBashforth<scalar_t, 3> AB3() {
    Eigen::Matrix<scalar_t, 3, 1> b(3); b << 5./12., -4./3., 23./12.;
    return {b};
}
 
template <class scalar_t = double>
static AdamsBashforth<scalar_t, 4> AB4() {
    Eigen::Matrix<scalar_t, 4, 1> b(4); b << -3./8., 37./24., -59./24., 55./24.;
    return {b};
}
 
template <class scalar_t = double>
static AdamsBashforth<scalar_t, 5> AB5() {
    Eigen::Matrix<scalar_t, 5, 1> b(5);
    b << 251./720., -637./360., 109./30., -1387./360., 1901./720.;
    return {b};
}
 
// Internal namespace containing implementation details for Explicit::of_order function.
namespace of_order_internal {
template <int order, class scalar_t>
struct _of_order_impl {
    static AdamsBashforth<scalar_t, order> of_order();
};
}  // namespace of_order_internal
 
template <int order, class scalar_t = double>
static AdamsBashforth<scalar_t, order> of_order() {
    return of_order_internal::_of_order_impl<order, scalar_t>::of_order();
}
 
}  // namespace ExplicitMultistep
}  // namespace integrators
}  // namespace mm
 
#endif  // MEDUSA_BITS_INTEGRATORS_ADAMSBASHFORTH_FWD_HPP_