From 0f04e97fa748b6740da4c9512b596d7d3a2788c5 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Wed, 6 Oct 2010 16:47:04 -0400 Subject: Add the Euler method. --- src/vZ/Euler.hpp | 78 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/vZ/Integrator.hpp | 51 +++++++++++++++++---------------- src/vZ/RK.hpp | 66 +++++++++++++++++++++++++++++++++++++++---- src/vZ/Simple.hpp | 25 +++++++++-------- 4 files changed, 179 insertions(+), 41 deletions(-) create mode 100644 src/vZ/Euler.hpp (limited to 'src/vZ') 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 * + * * + * 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 * + * . * + *************************************************************************/ + +#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 + class GenericEulerIntegrator : public GenericSimpleIntegrator + { + public: + typedef typename GenericSimpleIntegrator::Scalar Scalar; + typedef typename GenericSimpleIntegrator::Function Function; + + GenericEulerIntegrator(Function f) + : GenericSimpleIntegrator(f, s_a, s_b) { } + ~GenericEulerIntegrator() { } + + private: + typedef typename GenericSimpleIntegrator::ACoefficients ACoefficients; + typedef typename GenericSimpleIntegrator::BCoefficients BCoefficients; + + static ACoefficients s_a; + static BCoefficients s_b; + + static Scalar s_bArr[1]; + }; + + // Type alias + typedef GenericEulerIntegrator EulerIntegrator; + + // Implementation + + template + typename GenericEulerIntegrator::Scalar + GenericEulerIntegrator::s_bArr[1] = { + Scalar(1) + }; + + template + typename GenericEulerIntegrator::ACoefficients + GenericEulerIntegrator::s_a; + + template + typename GenericEulerIntegrator::BCoefficients + GenericEulerIntegrator::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 + // y' = f(x, y); y(x0) = y0 + // then an Integrator could be constructed as Integrator(f, dt).y(y0).x(x0) + template class GenericIntegrator { public: - typedef std::tr1::function Function; + typedef typename Traits::Scalar Scalar; + typedef std::tr1::function 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 + template void - GenericIntegrator::integrate(T t_final) + GenericIntegrator::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 - class GenericRKIntegrator : public GenericIntegrator + template + class GenericRKIntegrator : public GenericIntegrator { public: - typedef typename GenericIntegrator::Function Function; + typedef typename GenericIntegrator::Scalar Scalar; + typedef typename GenericIntegrator::Function Function; protected: // Coefficients in the tableau representation of the RK algorithm - typedef std::vector > ACoefficients; - typedef std::vector BCoefficients; + typedef std::vector > ACoefficients; + typedef std::vector BCoefficients; - GenericRKIntegrator(Function f, T dt) : GenericIntegrator(f, dt) { } + // Result vectors + typedef std::vector KVector; + + GenericRKIntegrator(Function f) : GenericIntegrator(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 RKIntegrator; + + // Implementation + + template + typename GenericRKIntegrator::KVector + GenericRKIntegrator::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::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 + Y + GenericRKIntegrator::calculateY(const KVector& k, const BCoefficients& b) + const + { + Y y = this->y(); + + for (typename std::vector::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 - class GenericSimpleIntegrator : public GenericRKIntegrator + template + class GenericSimpleIntegrator : public GenericRKIntegrator { public: - typedef typename GenericIntegrator::Function Function; + typedef typename GenericRKIntegrator::Scalar Scalar; + typedef typename GenericRKIntegrator::Function Function; protected: - typedef typename GenericRKIntegrator::ACoefficients ACoefficients; - typedef typename GenericRKIntegrator::BCoefficients BCoefficients; + typedef typename GenericRKIntegrator::ACoefficients ACoefficients; + typedef typename GenericRKIntegrator::BCoefficients BCoefficients; + typedef typename GenericRKIntegrator::KVector KVector; - GenericSimpleIntegrator(Function f, T dt, - ACoefficients a, BCoefficients b) - : GenericIntegrator(f, dt), m_a(a), m_b(b) { } + GenericSimpleIntegrator(Function f, ACoefficients a, BCoefficients b) + : GenericRKIntegrator(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 + template void - GenericSimpleIntegrator::step(T& t, T& dt) + GenericSimpleIntegrator::step() { + this->y(calculateY(calculateK(m_a), m_b)); + this->x(this->x() + this->h()); } } -- cgit v1.2.3