From 2829763fc5f9d222a966402bf073083dbc1da51c Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Wed, 13 Oct 2010 15:45:40 -0400
Subject: Implement FSAL.

---
 src/vZ/RK.hpp | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

(limited to 'src/vZ/RK.hpp')

diff --git a/src/vZ/RK.hpp b/src/vZ/RK.hpp
index 97ab487..8964741 100644
--- a/src/vZ/RK.hpp
+++ b/src/vZ/RK.hpp
@@ -46,6 +46,8 @@ namespace vZ
 
     // Perform the stages of an RK integration
     KVector calculateK(const ACoefficients& a) const;
+    KVector calculateK(Y& y, const ACoefficients& a) const;
+    KVector calculateK(Y k1, Y& y, const ACoefficients& a) const;
     Y calculateY(const KVector& k, const BCoefficients& b) const;
   };
 
@@ -57,12 +59,27 @@ namespace vZ
   template <typename Y>
   typename GenericRKIntegrator<Y>::KVector
   GenericRKIntegrator<Y>::calculateK(const ACoefficients& a) const
+  {
+    Y ytemp; // Ignored
+    return calculateK(ytemp, a);
+  }
+
+  template <typename Y>
+  typename GenericRKIntegrator<Y>::KVector
+  GenericRKIntegrator<Y>::calculateK(Y& y, const ACoefficients& a) const
+  {
+    return calculateK(this->f()(this->x(), this->y()), y, a);
+  }
+
+  template <typename Y>
+  typename GenericRKIntegrator<Y>::KVector
+  GenericRKIntegrator<Y>::calculateK(Y k1, Y& y, const ACoefficients& a) const
   {
     KVector k;
     k.reserve(a.size() + 1);
 
     // k1
-    k.push_back(this->f()(this->x(), this->y()));
+    k.push_back(k1);
 
     // k2..n
     for (typename ACoefficients::const_iterator i = a.begin();
@@ -70,7 +87,7 @@ namespace vZ
          ++i)
     {
       Scalar c(0);
-      Y y = this->y();
+      y = this->y();
       for (typename std::vector<Scalar>::size_type j = 0; j < i->size(); ++j) {
         Scalar aij = i->at(j);
         c += aij;
-- 
cgit v1.2.3