RKExplicit_fwd.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_INTEGRATORS_RKEXPLICIT_FWD_HPP_
#define MEDUSA_BITS_INTEGRATORS_RKEXPLICIT_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;
 
template <typename Scalar, int num_stages>
class RKExplicit {
  public:
    typedef Scalar scalar_t;  
    enum {  stages = num_stages };
 
  private:
    Eigen::Matrix<scalar_t, stages, 1> alpha_;  
    Eigen::Matrix<scalar_t, stages, stages> beta_;  
    Eigen::Matrix<scalar_t, stages, 1> gamma_;  
 
  public:
    RKExplicit(const Eigen::Matrix<scalar_t, stages, 1>& alpha,
               const Eigen::Matrix<scalar_t, stages, stages>& beta,
               const Eigen::Matrix<scalar_t, stages, 1>& gamma)
            : alpha_(alpha), beta_(beta), gamma_(gamma) {}
 
    template <typename func_t>
    class Integrator {
        const RKExplicit<scalar_t, num_stages> method_;  
        scalar_t t0_;  
        scalar_t dt_;  
        int max_steps;  
        Eigen::VectorXd y0_;  
        const func_t& func_;  
 
      public:
        Integrator(const RKExplicit<scalar_t, num_stages>& method, const func_t& func, scalar_t t0,
                   scalar_t t_max, scalar_t dt, Eigen::VectorXd y0) :
                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::VectorXd y;  
            int cur_step;  
 
          public:
            explicit IterationStep(const Integrator& integrator)
                    : integrator(integrator), t(integrator.t0_), y(integrator.y0_), cur_step(0) {}
 
            IterationStep(const Integrator& integrator, int) : integrator(integrator), t(), y(),
                                                               cur_step(integrator.max_steps + 1) {}
 
            IterationStep& operator++() {
                y = integrator.method_.step(integrator.func_, t, y, integrator.dt_);
                t += integrator.dt_;
                cur_step++;
                return *this;
            }
 
            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; }
 
            Eigen::VectorXd value() const { return y; }
 
            scalar_t& time() { return t; }
 
            Eigen::VectorXd& value() { return y; }
 
            friend std::ostream& operator<<(std::ostream& os, const IterationStep& step) {
                return os << "t = " << step.t << "; y = " << step.y;
            }
        };
 
        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 << "Explicit RK 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>(*this, func, t0, t_max, dt, y0);
    }
 
  private:
    template <typename func_t>
    Eigen::VectorXd step(const func_t& func, scalar_t t, const Eigen::VectorXd& y,
                         scalar_t dt) const {
        Eigen::Matrix<scalar_t, Eigen::Dynamic, stages> k(y.size(), static_cast<int>(stages));
        k.col(0) = func(t + alpha_(0) * dt, y);
        for (int i = 1; i < stages; ++i) {
            k.col(i) = func(t + alpha_(i) * dt,
                            y + dt * k.leftCols(i) * beta_.row(i).head(i).transpose());
        }
        return y + dt * k * gamma_;
    }
 
  public:
    template <typename scalar_t, int stages>
    friend std::ostream& operator<<(std::ostream& os, const RKExplicit<scalar_t, stages>&);
 
    template <typename Scalar2, int num_steps2>
    friend class AdamsBashforth;
};
 
template <typename scalar_t, int num_stages>
std::ostream& operator<<(std::ostream& os, const RKExplicit<scalar_t, num_stages>& method) {
    os << "RKExplicit with " << num_stages << " stages and "
       << "alpha = " << method.alpha_ << ", "
       << "beta = " << method.beta_ << ", "
       << "gamma = " << method.gamma_;
    return os;
}
 
namespace integrators {
 
namespace Explicit {
 
template <class scalar_t = double>
RKExplicit<scalar_t, 4> RK4();
 
template <class scalar_t = double>
RKExplicit<scalar_t, 4> RK38();
 
template <class scalar_t = double>
RKExplicit<scalar_t, 1> Euler();
 
template <class scalar_t = double>
RKExplicit<scalar_t, 2> Midpoint();
 
template <class scalar_t = double>
RKExplicit<scalar_t, 3> RK3();
 
template <class scalar_t = double>
RKExplicit<scalar_t, 6> Fehlberg5();
 
namespace of_order_internal {
template <int order, class scalar_t>
struct _of_order_impl { static RKExplicit<scalar_t, order> of_order(); };
}
 
template <int order, class scalar_t = double>
static RKExplicit<scalar_t, order> of_order() {
    return of_order_internal::_of_order_impl<order, scalar_t>::of_order();
}
 
}  // namespace Explicit
}  // namespace integrators
}  // namespace mm
 
#endif  // MEDUSA_BITS_INTEGRATORS_RKEXPLICIT_FWD_HPP_