From f123be63a5760782c563d81d7842de6cc3e5646b Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Mon, 11 Oct 2010 02:51:41 -0400
Subject: Support adaptive integration with equation systems.

---
 src/vZ/EquationSystem.hpp | 14 ++++++++++++++
 tests/EquationSystem.cpp  | 10 ++++++----
 2 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/vZ/EquationSystem.hpp b/src/vZ/EquationSystem.hpp
index 0963baf..6409f18 100644
--- a/src/vZ/EquationSystem.hpp
+++ b/src/vZ/EquationSystem.hpp
@@ -21,6 +21,8 @@
 #ifndef VZ_EQUATIONSYSTEM_HPP
 #define VZ_EQUATIONSYSTEM_HPP
 
+#include <algorithm>
+#include <cmath>
 #include <cstddef>
 
 namespace vZ
@@ -149,6 +151,18 @@ namespace vZ
     }
     return *this;
   }
+
+  template <std::size_t N, typename T>
+  typename EquationSystem<N, T>::Scalar
+  abs(const EquationSystem<N, T>& es)
+  {
+    typename EquationSystem<N, T>::Scalar ret(0);
+    for (std::size_t i = 0; i < N; ++i) {
+      using std::abs;
+      ret = std::max(ret, abs(es[i]));
+    }
+    return ret;
+  }
 }
 
 #endif // VZ_EQUATIONSYSTEM_HPP
diff --git a/tests/EquationSystem.cpp b/tests/EquationSystem.cpp
index 973dd30..497a065 100644
--- a/tests/EquationSystem.cpp
+++ b/tests/EquationSystem.cpp
@@ -22,8 +22,8 @@ main()
   Y y;
   y[0] = 1.0;
   y[1] = 1.0;
-  vZ::GenericEulerIntegrator<Y> integrator(f);
-  integrator.y(y).x(0.0).h(0.01);
+  vZ::GenericDP45Integrator<Y> integrator(f);
+  integrator.tol(1e-6).y(y).x(0.0).h(0.06);
 
   integrator.integrate(2.0);
 
@@ -33,10 +33,12 @@ main()
   std::cout << std::setprecision(10)
             << "Numerical:  " << actual << std::endl
             << "Expected:   " << expected  << std::endl
-            << "Iterations: " << integrator.iterations() << std::endl;
+            << "h:          " << integrator.h() << std::endl
+            << "Iterations: " << integrator.iterations() << std::endl
+            << "Rejections: " << integrator.rejections() << std::endl;
 
   double error = std::fabs(expected - actual)/expected;
-  if (error > 0.01) {
+  if (error > 6.0e-7) {
     std::cerr << "Error:      " << 100.0*error << "%" << std::endl;
     return EXIT_FAILURE;
   } else {
-- 
cgit v1.2.3