euler.hh

#ifndef VTESTBED_CORE_EULER_HH
#define VTESTBED_CORE_EULER_HH

#include <algorithm>

#include <vtestbed/base/blitz.hh>

namespace vtb {

    namespace core {

        template <class T>
        class Euler_solver {

        public:
            using value_type = T;

        private:
            T _minstep = T(1e-3);
            T _step{1e-3};
            bool _verbose = false;

        public:

            Euler_solver() = default;
            Euler_solver(const Euler_solver&) = default;
            Euler_solver& operator=(const Euler_solver&) = default;
            Euler_solver(Euler_solver&&) = default;
            Euler_solver& operator=(Euler_solver&&) = default;
            ~Euler_solver() = default;

            template <class Function, int N>
            inline blitz::TinyVector<T,N>
            solve(Function f, T t0, T t1, const blitz::TinyVector<T,N>& x0) {
                using vec = blitz::TinyVector<T,N>;
                using std::min;
                bool finished = false;
                const T h_min = this->_minstep;
                T h = this->_step;
                T t_n = t0;
                vec x_n = x0;
                while (!finished) {
                    if (this->_verbose) {
                        std::clog
                            << std::setw(8) << "euler"
                            << std::setw(12) << t_n
                            << std::setw(12) << h
                            << std::setw(10) << x_n
                            << std::endl;
                    }
                    if (t_n+h+h_min > t1) {
                        h = t1-t_n;
                        finished = true;
                    }
                    x_n += h*f(t_n, x_n);
                    t_n += h;
                }
                if (this->_verbose) {
                    std::clog
                        << std::setw(8) << "euler"
                        << std::setw(12) << t_n
                        << std::setw(12) << h
                        << std::setw(10) << x_n
                        << std::endl;
                }
                return x_n;
            }

            template <class Function, int N>
            inline blitz::TinyVector<T,N>
            operator()(Function f, T t0, T t1, const blitz::TinyVector<T,N>& x0) {
                return this->solve<Function,N>(f, t0, t1, x0);
            }

            inline void
            step(T rhs) {
                if (!(rhs > T(0))) {
                    throw std::invalid_argument("Euler_solver::step <= 0");
                }
                this->_step = rhs;
            }

            inline void
            min_step(T rhs) {
                if (!(rhs > T(0))) {
                    throw std::invalid_argument("Euler_solver::min_step <= 0");
                }
                this->_minstep = rhs;
            }

            inline T min_step() const { return this->_minstep; }
            inline T step() const { return this->_step; }
            inline void verbose(bool rhs) { this->_verbose = rhs; }
            inline bool verbose() const { return this->_verbose; }

        };


    }

}

#endif // vim:filetype=cpp