Changed rk4 signature.

matrix/integration.hpp

@@ -7,19 +7,19 @@ namespace matrix {
 template<typename Type, size_t M>
 int integrate_rk4(
     Vector<Type, M> (*f)(Type, Vector<Type, M>),
-    Vector<Type, M> & y,
-    Type & t,
-    Type h
+    const Vector<Type, M> & y0,
+    Type t0,
+    Type h,
+    Vector<Type, M> & y1
 )
 {
     // https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
     Vector<Type, M> k1, k2, k3, k4;
-    k1 = f(t, y);
-    k2 = f(t + h/2, y + k1*h/2);
-    k3 = f(t + h/2, y + k2*h/2);
-    k4 = f(t + h, y + k3*h);
-    y += (k1 + k2*2 + k3*2 + k4)*(h/6);
-    t += h;
+    k1 = f(t0, y0);
+    k2 = f(t0 + h/2, y0 + k1*h/2);
+    k3 = f(t0 + h/2, y0 + k2*h/2);
+    k4 = f(t0 + h, y0 + k3*h);
+    y1 = y0 + (k1 + k2*2 + k3*2 + k4)*(h/6);
     return 0;
 }
 

test/integration.cpp

@@ -8,9 +8,10 @@ using namespace matrix;
 // instantiate template to ensure coverage check
 template int integrate_rk4<float, 6>(
     Vector<float, 6> (*f)(float, Vector<float, 6>),
-    Vector<float, 6> & y,
-    float & t,
-    float h
+    const Vector<float, 6> & y0,
+    float t0,
+    float h,
+    Vector<float, 6> & y1
 );
 
 Vector<float, 6> f(float t, Vector<float, 6> y);
@@ -22,12 +23,12 @@ Vector<float, 6> f(float t, Vector<float, 6> y) {
 int main()
 {
     Vector<float, 6> y = ones<float, 6, 1>();
-    float h = 1;
     float t = 1;
+    float h = 0.1f;
     y.T().print();
-    integrate_rk4(f, y, t, h);
+    integrate_rk4(f, y, t, h, y);
     y.T().print();
-    assert(y == (ones<float, 6, 1>()*2));
+    assert(y == (ones<float, 6, 1>()*1.1f));
     return 0;
 }