diff options
Diffstat (limited to 'src/vZ')
| -rw-r--r-- | src/vZ/Euler.hpp | 78 | ||||
| -rw-r--r-- | src/vZ/Integrator.hpp | 51 | ||||
| -rw-r--r-- | src/vZ/RK.hpp | 66 | ||||
| -rw-r--r-- | src/vZ/Simple.hpp | 25 |
4 files changed, 179 insertions, 41 deletions
diff --git a/src/vZ/Euler.hpp b/src/vZ/Euler.hpp new file mode 100644 index 0000000..56c5678 --- /dev/null +++ b/src/vZ/Euler.hpp @@ -0,0 +1,78 @@ +/************************************************************************* + * Copyright (C) 2009-2010 Tavian Barnes <tavianator@gmail.com> * + * * + * This file is part of The vZ Library. * + * * + * The vZ Library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Lesser General Public License as * + * published by the Free Software Foundation; either version 3 of the * + * License, or (at your option) any later version. * + * * + * The vZ Library is distributed in the hope that it will be useful, but * + * WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * Lesser General Public License for more details. * + * * + * You should have received a copy of the GNU Lesser General Public * + * License along with this program. If not, see * + * <http://www.gnu.org/licenses/>. * + *************************************************************************/ + +#ifndef VZ_EULER_HPP +#define VZ_EULER_HPP + +namespace vZ +{ + // Euler method + // + // Simplest Runge-Kutta method + // First order + // Its tableau is: + // + // 0| + // -+- + // |1 + // + // y[n + 1] = y[n] + dt*f(y[n]) + template <typename Y> + class GenericEulerIntegrator : public GenericSimpleIntegrator<Y> + { + public: + typedef typename GenericSimpleIntegrator<Y>::Scalar Scalar; + typedef typename GenericSimpleIntegrator<Y>::Function Function; + + GenericEulerIntegrator(Function f) + : GenericSimpleIntegrator<Y>(f, s_a, s_b) { } + ~GenericEulerIntegrator() { } + + private: + typedef typename GenericSimpleIntegrator<Y>::ACoefficients ACoefficients; + typedef typename GenericSimpleIntegrator<Y>::BCoefficients BCoefficients; + + static ACoefficients s_a; + static BCoefficients s_b; + + static Scalar s_bArr[1]; + }; + + // Type alias + typedef GenericEulerIntegrator<double> EulerIntegrator; + + // Implementation + + template <typename Y> + typename GenericEulerIntegrator<Y>::Scalar + GenericEulerIntegrator<Y>::s_bArr[1] = { + Scalar(1) + }; + + template <typename Y> + typename GenericEulerIntegrator<Y>::ACoefficients + GenericEulerIntegrator<Y>::s_a; + + template <typename Y> + typename GenericEulerIntegrator<Y>::BCoefficients + GenericEulerIntegrator<Y>::s_b(s_bArr, s_bArr + 1); +} + +#endif // VZ_EULER_HPP diff --git a/src/vZ/Integrator.hpp b/src/vZ/Integrator.hpp index 25a84f6..539c313 100644 --- a/src/vZ/Integrator.hpp +++ b/src/vZ/Integrator.hpp @@ -30,37 +30,40 @@ namespace vZ // // All integration methods derrive from this class // If the initial value problem is specified as - // y' = f(t, y); y(t0) = y0 - // then an Integrator could be constructed as Integrator(f, dt).y(y0).t(t0) - template <typename T> + // y' = f(x, y); y(x0) = y0 + // then an Integrator could be constructed as Integrator(f, dt).y(y0).x(x0) + template <typename Y> class GenericIntegrator { public: - typedef std::tr1::function<T (T, T)> Function; + typedef typename Traits<Y>::Scalar Scalar; + typedef std::tr1::function<Y (Scalar, Y)> Function; - // By default, y and t start at zero - GenericIntegrator(Function f, T dt) - : m_f(f), m_y(0), m_t(0), m_dt(dt) { } + // By default, y and t start at zero, h starts UNDEFINED + GenericIntegrator(Function f) + : m_f(f), m_y(0), m_x(0), m_h() { } virtual ~GenericIntegrator() { } - GenericIntegrator& y(T y) { m_y = y; return *this; } - GenericIntegrator& t(T t) { m_t = t; return *this; } - GenericIntegrator& dt(T dt) { m_dt = dt; return *this; } + GenericIntegrator& y(Y y) { m_y = y; return *this; } + GenericIntegrator& x(Scalar x) { m_x = x; return *this; } + GenericIntegrator& h(Scalar h) { m_h = h; return *this; } - T y() const { return m_y; } - T t() const { return m_t; } - T dt() const { return m_dt; } + Y y() const { return m_y; } + Scalar x() const { return m_x; } + Scalar h() const { return m_h; } - // Integrate until time t - void integrate(T t_final); + // Integrate until x == x_final + void integrate(Scalar x_final); protected: - virtual void step(T& t, T& dt) = 0; - Function m_f; + virtual void step() = 0; + + const Function& f() const { return m_f; } private: - T m_y; - T m_t, m_dt; + Function m_f; + Y m_y; + Scalar m_x, m_h; }; // Type alias @@ -68,13 +71,13 @@ namespace vZ // Implementations - template <typename T> + template <typename Y> void - GenericIntegrator<T>::integrate(T t_final) + GenericIntegrator<Y>::integrate(Scalar x_final) { - while (m_t < t_final) { - m_dt = std::min(m_dt, t_final - m_t); - step(m_t, m_dt); + while (m_x < x_final) { + m_h = std::min(m_h, x_final - m_x); + step(); } } } diff --git a/src/vZ/RK.hpp b/src/vZ/RK.hpp index 3ba4b8a..df608ff 100644 --- a/src/vZ/RK.hpp +++ b/src/vZ/RK.hpp @@ -26,23 +26,77 @@ namespace vZ { // Base class for Runge-Kutta type algorithms - template <typename T> - class GenericRKIntegrator : public GenericIntegrator<T> + template <typename Y> + class GenericRKIntegrator : public GenericIntegrator<Y> { public: - typedef typename GenericIntegrator<T>::Function Function; + typedef typename GenericIntegrator<Y>::Scalar Scalar; + typedef typename GenericIntegrator<Y>::Function Function; protected: // Coefficients in the tableau representation of the RK algorithm - typedef std::vector<std::vector<T> > ACoefficients; - typedef std::vector<T> BCoefficients; + typedef std::vector<std::vector<Scalar> > ACoefficients; + typedef std::vector<Scalar> BCoefficients; - GenericRKIntegrator(Function f, T dt) : GenericIntegrator<T>(f, dt) { } + // Result vectors + typedef std::vector<Y> KVector; + + GenericRKIntegrator(Function f) : GenericIntegrator<Y>(f) { } virtual ~GenericRKIntegrator() { } + + // Perform the stages of an RK integration + KVector calculateK(const ACoefficients& a) const; + Y calculateY(const KVector& k, const BCoefficients& b) const; }; // Type alias typedef GenericRKIntegrator<double> RKIntegrator; + + // Implementation + + template <typename Y> + typename GenericRKIntegrator<Y>::KVector + GenericRKIntegrator<Y>::calculateK(const ACoefficients& a) const + { + KVector k; + k.reserve(a.size() + 1); + + // k1 + k.push_back(this->f()(this->x(), this->y())); + + // k2..n + for (typename ACoefficients::const_iterator i = a.begin(); + i != a.end(); + ++i) + { + Scalar c = 0; + Y y = this->y(); + for (typename std::vector<Scalar>::size_type j = 0; j < i->size(); ++j) { + Scalar aij = i->at(j); + c += aij; + + y += aij*this->h()*k.at(j); + } + + k.push_back(this->f()(this->x() + c, y)); + } + + return k; + } + + template <typename Y> + Y + GenericRKIntegrator<Y>::calculateY(const KVector& k, const BCoefficients& b) + const + { + Y y = this->y(); + + for (typename std::vector<Scalar>::size_type i = 0; i < k.size(); ++i) { + y += this->h()*b.at(i)*k.at(i); + } + + return y; + } } #endif // VZ_RK_HPP diff --git a/src/vZ/Simple.hpp b/src/vZ/Simple.hpp index 103357e..9d9be46 100644 --- a/src/vZ/Simple.hpp +++ b/src/vZ/Simple.hpp @@ -26,22 +26,23 @@ namespace vZ { // Base class for non-adaptive RK-style algorithms - template <typename T> - class GenericSimpleIntegrator : public GenericRKIntegrator<T> + template <typename Y> + class GenericSimpleIntegrator : public GenericRKIntegrator<Y> { public: - typedef typename GenericIntegrator<T>::Function Function; + typedef typename GenericRKIntegrator<Y>::Scalar Scalar; + typedef typename GenericRKIntegrator<Y>::Function Function; protected: - typedef typename GenericRKIntegrator<T>::ACoefficients ACoefficients; - typedef typename GenericRKIntegrator<T>::BCoefficients BCoefficients; + typedef typename GenericRKIntegrator<Y>::ACoefficients ACoefficients; + typedef typename GenericRKIntegrator<Y>::BCoefficients BCoefficients; + typedef typename GenericRKIntegrator<Y>::KVector KVector; - GenericSimpleIntegrator(Function f, T dt, - ACoefficients a, BCoefficients b) - : GenericIntegrator<T>(f, dt), m_a(a), m_b(b) { } + GenericSimpleIntegrator(Function f, ACoefficients a, BCoefficients b) + : GenericRKIntegrator<Y>(f), m_a(a), m_b(b) { } virtual ~GenericSimpleIntegrator() { } - virtual void step(T& t, T& dt); + void step(); private: ACoefficients m_a; @@ -53,10 +54,12 @@ namespace vZ // Implementations - template <typename T> + template <typename Y> void - GenericSimpleIntegrator<T>::step(T& t, T& dt) + GenericSimpleIntegrator<Y>::step() { + this->y(calculateY(calculateK(m_a), m_b)); + this->x(this->x() + this->h()); } } |