summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorTavian Barnes <tavianator@gmail.com>2010-10-06 16:47:04 -0400
committerTavian Barnes <tavianator@gmail.com>2010-10-06 16:47:04 -0400
commit0f04e97fa748b6740da4c9512b596d7d3a2788c5 (patch)
treeeb03f2de5798c78e13c53b305758bf727403f872
parent2e25da27f14566000fb34d3859bfb470bf5fd1da (diff)
downloadvz-0f04e97fa748b6740da4c9512b596d7d3a2788c5.tar.xz
Add the Euler method.
-rw-r--r--src/Makefile.am1
-rw-r--r--src/vZ.hpp1
-rw-r--r--src/vZ/Euler.hpp78
-rw-r--r--src/vZ/Integrator.hpp51
-rw-r--r--src/vZ/RK.hpp66
-rw-r--r--src/vZ/Simple.hpp25
-rw-r--r--tests/Euler.cpp25
-rw-r--r--tests/Makefile.am4
-rw-r--r--tests/compile.cpp8
9 files changed, 208 insertions, 51 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 5ea10b7..76c30b1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -18,6 +18,7 @@
###########################################################################
nobase_include_HEADERS = vZ.hpp \
+ vZ/Euler.hpp \
vZ/Integrator.hpp \
vZ/RK.hpp \
vZ/Simple.hpp \
diff --git a/src/vZ.hpp b/src/vZ.hpp
index 64c894d..12e4cbc 100644
--- a/src/vZ.hpp
+++ b/src/vZ.hpp
@@ -27,5 +27,6 @@
#include <vZ/Integrator.hpp>
#include <vZ/RK.hpp>
#include <vZ/Simple.hpp>
+#include <vZ/Euler.hpp>
#endif // VZ_HPP
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());
}
}
diff --git a/tests/Euler.cpp b/tests/Euler.cpp
new file mode 100644
index 0000000..5666cd8
--- /dev/null
+++ b/tests/Euler.cpp
@@ -0,0 +1,25 @@
+#include "vZ.hpp"
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+
+// y' = y (y == C*exp(t))
+double
+f(double t, double y)
+{
+ return y;
+}
+
+int
+main()
+{
+ vZ::EulerIntegrator integrator(f);
+ integrator.y(1.0).x(0.0).h(0.01);
+
+ integrator.integrate(2.0);
+
+ std::cout << integrator.y() << std::endl
+ << std::exp(2.0) << std::endl;
+
+ return EXIT_SUCCESS;
+}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 156bc06..a3d5a77 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -19,7 +19,7 @@
INCLUDES = -I$(top_srcdir)/src
-check_PROGRAMS = compile-test
+check_PROGRAMS = Euler-test
TESTS = $(check_PROGRAMS)
-compile_test_SOURCES = compile.cpp
+Euler_test_SOURCES = Euler.cpp
diff --git a/tests/compile.cpp b/tests/compile.cpp
deleted file mode 100644
index 17381b4..0000000
--- a/tests/compile.cpp
+++ /dev/null
@@ -1,8 +0,0 @@
-#include "vZ.hpp"
-#include <cstdlib>
-
-int
-main()
-{
- return EXIT_SUCCESS;
-}