Add the Euler method.

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());
   }
 }