From 0f04e97fa748b6740da4c9512b596d7d3a2788c5 Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Wed, 6 Oct 2010 16:47:04 -0400
Subject: Add the Euler method.

---
 src/vZ/Integrator.hpp | 51 +++++++++++++++++++++++++++------------------------
 1 file changed, 27 insertions(+), 24 deletions(-)

(limited to 'src/vZ/Integrator.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();
     }
   }
 }
-- 
cgit v1.2.3