Add Heun-Euler adaptive method.

11 files changed, 276 insertions, 17 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 9ad815b..04dabb8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -18,7 +18,9 @@
 ###########################################################################
 
 nobase_include_HEADERS = vZ.hpp                                                \
+                         vZ/Adaptive.hpp                                       \
                          vZ/Euler.hpp                                          \
+                         vZ/HE12.hpp                                           \
                          vZ/Heun.hpp                                           \
                          vZ/Integrator.hpp                                     \
                          vZ/Midpoint.hpp                                       \
diff --git a/src/vZ.hpp b/src/vZ.hpp
index 7a311b1..9ce7beb 100644
--- a/src/vZ.hpp
+++ b/src/vZ.hpp
@@ -31,5 +31,7 @@
 #include <vZ/Midpoint.hpp>
 #include <vZ/Heun.hpp>
 #include <vZ/RK4.hpp>
+#include <vZ/Adaptive.hpp>
+#include <vZ/HE12.hpp>
 
 #endif // VZ_HPP
diff --git a/src/vZ/Adaptive.hpp b/src/vZ/Adaptive.hpp
new file mode 100644
index 0000000..864b038
--- /dev/null
+++ b/src/vZ/Adaptive.hpp
@@ -0,0 +1,109 @@
+/*************************************************************************
+ * 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_ADAPTIVE_HPP
+#define VZ_ADAPTIVE_HPP
+
+#include <cmath>
+
+namespace vZ
+{
+  // Base class for adaptive RK-style algorithms
+  template <typename Y>
+  class GenericAdaptiveIntegrator : public GenericRKIntegrator<Y>
+  {
+  public:
+    typedef typename GenericRKIntegrator<Y>::Scalar   Scalar;
+    typedef typename GenericRKIntegrator<Y>::Function Function;
+
+    GenericAdaptiveIntegrator& atol(Scalar tol) { m_atol = tol; return *this; }
+    GenericAdaptiveIntegrator& rtol(Scalar tol) { m_rtol = tol; return *this; }
+
+    Scalar atol() const { return m_atol; }
+    Scalar rtol() const { return m_rtol; }
+
+  protected:
+    typedef typename GenericRKIntegrator<Y>::ACoefficients ACoefficients;
+    typedef typename GenericRKIntegrator<Y>::BCoefficients BCoefficients;
+    typedef typename GenericRKIntegrator<Y>::KVector       KVector;
+
+    GenericAdaptiveIntegrator(Function f, unsigned int order,
+                              ACoefficients a, BCoefficients b,
+                              BCoefficients bStar)
+      : GenericRKIntegrator<Y>(f), m_order(order),
+        m_a(a), m_b(b), m_bStar(bStar)
+    { }
+    virtual ~GenericAdaptiveIntegrator() { }
+
+    void step();
+
+  private:
+    Scalar m_atol, m_rtol;
+    unsigned int m_order;
+    ACoefficients m_a;
+    BCoefficients m_b, m_bStar;
+  };
+
+  // Type alias
+  typedef GenericAdaptiveIntegrator<double> AdaptiveIntegrator;
+
+  // Implementations
+
+  template <typename Y>
+  void
+  GenericAdaptiveIntegrator<Y>::step()
+  {
+    static const Scalar S = Scalar(19)/Scalar(20); // Arbitrary saftey factor
+    Scalar newH;
+    Y y;
+
+    // Attempt the integration step in a loop
+    bool rejected = true;
+    while (rejected) {
+      KVector k = calculateK(m_a);
+      y         = calculateY(k, m_b);
+      Y yStar   = calculateY(k, m_bStar);
+
+      // Get an error estimate
+      using std::abs;
+      using std::pow;
+      Scalar delta = abs(y - yStar);
+      Scalar scale = m_atol + std::max(abs(y), abs(this->y()))*m_rtol;
+
+      newH = S*this->h()*pow(scale/delta, Scalar(1)/m_order);
+
+      if (delta > scale) {
+        // Reject the step
+        this->h(newH);
+      } else {
+        rejected = false;
+      }
+    }
+
+    // Update x and y
+    this->y(y);
+    this->x(this->x() + this->h());
+
+    // Adjust the stepsize for the next iteration
+    this->h(newH);
+  }
+}
+
+#endif // VZ_ADAPTIVE_HPP
diff --git a/src/vZ/Euler.hpp b/src/vZ/Euler.hpp
index 56c5678..033da2e 100644
--- a/src/vZ/Euler.hpp
+++ b/src/vZ/Euler.hpp
@@ -29,9 +29,9 @@ namespace vZ
   // First order
   // Its tableau is:
   //
-  //   0|
-  //   -+-
-  //    |1
+  //   0 |
+  //   --+--
+  //     | 1
   //
   // y[n + 1] = y[n] + dt*f(y[n])
   template <typename Y>
diff --git a/src/vZ/HE12.hpp b/src/vZ/HE12.hpp
new file mode 100644
index 0000000..0f7a299
--- /dev/null
+++ b/src/vZ/HE12.hpp
@@ -0,0 +1,107 @@
+/*************************************************************************
+ * 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_HE12_HPP
+#define VZ_HE12_HPP
+
+namespace vZ
+{
+  // Heun-Euler method
+  //
+  // Second-order with embedded first-order
+  // Its tableau is:
+  //
+  //   0 |
+  //   1 | 1
+  //   --+---------
+  //     | 1/2  1/2
+  //     |  1    0
+  //
+  //    k1    = dt*f(y[n])
+  //    k2    = dt*f(y[n] + dt*k1)
+  // y[n + 1] = y[n] + 1/2*(k1 + k2)
+  template <typename Y>
+  class GenericHE12Integrator : public GenericAdaptiveIntegrator<Y>
+  {
+  public:
+    typedef typename GenericAdaptiveIntegrator<Y>::Scalar   Scalar;
+    typedef typename GenericAdaptiveIntegrator<Y>::Function Function;
+
+    GenericHE12Integrator(Function f)
+      : GenericAdaptiveIntegrator<Y>(f, 2, s_a, s_b, s_bStar) { }
+    ~GenericHE12Integrator() { }
+
+  private:
+    typedef typename GenericAdaptiveIntegrator<Y>::ACoefficients ACoefficients;
+    typedef typename GenericAdaptiveIntegrator<Y>::BCoefficients BCoefficients;
+
+    static ACoefficients s_a;
+    static BCoefficients s_b;
+    static BCoefficients s_bStar;
+
+    static Scalar s_a2Arr[1];
+    static std::vector<Scalar> s_aArr[1];
+    static Scalar s_bArr[2];
+    static Scalar s_bStarArr[2];
+  };
+
+  // Type alias
+  typedef GenericHE12Integrator<double> HE12Integrator;
+
+  // Implementation
+
+  template <typename Y>
+  typename GenericHE12Integrator<Y>::Scalar
+  GenericHE12Integrator<Y>::s_a2Arr[1] = {
+    Scalar(1)
+  };
+
+  template <typename Y>
+  std::vector<typename GenericHE12Integrator<Y>::Scalar>
+  GenericHE12Integrator<Y>::s_aArr[1] = {
+    std::vector<Scalar>(s_a2Arr, s_a2Arr + 1)
+  };
+
+  template <typename Y>
+  typename GenericHE12Integrator<Y>::ACoefficients
+  GenericHE12Integrator<Y>::s_a(s_aArr, s_aArr + 1);
+
+  template <typename Y>
+  typename GenericHE12Integrator<Y>::Scalar
+  GenericHE12Integrator<Y>::s_bArr[2] = {
+    Scalar(1)/Scalar(2), Scalar(1)/Scalar(2)
+  };
+
+  template <typename Y>
+  typename GenericHE12Integrator<Y>::BCoefficients
+  GenericHE12Integrator<Y>::s_b(s_bArr, s_bArr + 2);
+
+  template <typename Y>
+  typename GenericHE12Integrator<Y>::Scalar
+  GenericHE12Integrator<Y>::s_bStarArr[2] = {
+    Scalar(1), Scalar(0)
+  };
+
+  template <typename Y>
+  typename GenericHE12Integrator<Y>::BCoefficients
+  GenericHE12Integrator<Y>::s_bStar(s_bStarArr, s_bStarArr + 2);
+}
+
+#endif // VZ_HE12_HPP
diff --git a/src/vZ/Heun.hpp b/src/vZ/Heun.hpp
index b70df9e..5e697af 100644
--- a/src/vZ/Heun.hpp
+++ b/src/vZ/Heun.hpp
@@ -32,7 +32,7 @@ namespace vZ
   //   0 |
   //   1 | 1
   //   --+---------
-  //     | 1/2  1/2
+  //     | 1/2 1/2
   //
   //    k1    = dt*f(y[n])
   //    k2    = dt*f(y[n] + dt*k1)
diff --git a/src/vZ/Midpoint.hpp b/src/vZ/Midpoint.hpp
index 2329457..a4bb906 100644
--- a/src/vZ/Midpoint.hpp
+++ b/src/vZ/Midpoint.hpp
@@ -29,10 +29,10 @@ namespace vZ
   // Two function evaluations per step
   // Its tableau is:
   //
-  //   0  |
-  //   1/2|1/2
-  //   ---+-----
-  //      |0   1
+  //   0   |
+  //   1/2 | 1/2
+  //   ----+------
+  //       | 0   1
   //
   //    k1    = dt*f(y[n])
   //    k2    = dt*f(y[n] + (dt/2)*k1)
diff --git a/src/vZ/RK4.hpp b/src/vZ/RK4.hpp
index fa50461..39fa306 100644
--- a/src/vZ/RK4.hpp
+++ b/src/vZ/RK4.hpp
@@ -29,12 +29,12 @@ namespace vZ
   // Four function evaluations per step
   // Its tableau is:
   //
-  //   0  |
-  //   1/2|1/2
-  //   1/2|0   1/2
-  //   1  |0   0   1
-  //   ---+---------------
-  //      |1/6 1/3 1/3 1/6
+  //   0   |
+  //   1/2 | 1/2
+  //   1/2 | 0   1/2
+  //   1   | 0   0   1
+  //   ----+----------------
+  //       | 1/6 1/3 1/3 1/6
   //
   //    k1    = dt*f(y[n])
   //    k2    = dt*f(y[n] + (dt/2)*k1)
diff --git a/src/vZ/Simple.hpp b/src/vZ/Simple.hpp
index 9d9be46..5e4f7f8 100644
--- a/src/vZ/Simple.hpp
+++ b/src/vZ/Simple.hpp
@@ -21,8 +21,6 @@
 #ifndef VZ_SIMPLE_HPP
 #define VZ_SIMPLE_HPP
 
-#include <vector>
-
 namespace vZ
 {
   // Base class for non-adaptive RK-style algorithms
diff --git a/tests/HE12.cpp b/tests/HE12.cpp
new file mode 100644
index 0000000..9a556ef
--- /dev/null
+++ b/tests/HE12.cpp
@@ -0,0 +1,39 @@
+#include "vZ.hpp"
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+
+// y' = y (y == C*exp(t))
+double
+f(double t, double y)
+{
+  return y;
+}
+
+int
+main()
+{
+  double tol = 1e-6;
+  vZ::HE12Integrator integrator(f);
+  integrator.atol(tol).rtol(tol).y(1.0).x(0.0).h(0.02);
+
+  integrator.integrate(2.0);
+
+  double actual   = integrator.y();
+  double expected = std::exp(2.0);
+
+  std::cout << std::setprecision(10)
+            << "Numerical: " << actual << std::endl
+            << "Expected:  " << expected  << std::endl
+            << "h:         " << integrator.h() << std::endl;
+
+  double error = std::fabs(expected - actual)/expected;
+  if (error > 8.7e-7) {
+    std::cerr << "Error:     " << 100.0*error << "%" << std::endl;
+    return EXIT_FAILURE;
+  } else {
+    std::cout << "Error:     " << 100.0*error << "%" << std::endl;
+    return EXIT_SUCCESS;
+  }
+}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 0baa463..e3c7601 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -22,10 +22,12 @@ INCLUDES = -I$(top_srcdir)/src
 check_PROGRAMS = Euler-test                                                    \
                  Midpoint-test                                                 \
                  Heun-test                                                     \
-                 RK4-test
+                 RK4-test                                                      \
+                 HE12-test
 TESTS          = $(check_PROGRAMS)
 
 Euler_test_SOURCES    = Euler.cpp
 Midpoint_test_SOURCES = Midpoint.cpp
 Heun_test_SOURCES     = Heun.cpp
 RK4_test_SOURCES      = RK4.cpp
+HE12_test_SOURCES     = HE12.cpp