Add an adaptive RK4 integrator

_code/pset_04/integrators/integrators.h

 #ifndef INTEGRATORS_H
 #define INTEGRATORS_H
 
+#include <utility>
+
 //--------------------------------------------------------------------------------------------------
 // A single step of Euler integration
 template <typename F, typename Vector, typename Scalar>
@@ -28,4 +30,32 @@ Vector rk4_step(F const& f, Vector const& position, Scalar step_size) {
     return result;
 }
 
+//--------------------------------------------------------------------------------------------------
+// A fourth order Runge-Kutta step that also returns a recommneded next step size, based on a local
+// error tolerance.
+template <typename F, typename Vector, typename Scalar>
+std::pair<Vector, Scalar> adaptive_rk4_step(F const& f, Vector const& position, Scalar step_size, Scalar tolerance) {
+    // It's key that bigger * smaller != 1, so that it's possible to reach any step size.
+    static constexpr Scalar bigger = 1.2;
+    static constexpr Scalar smaller = 0.8;
+
+    Vector const full_step = rk4_step(f, position, step_size);
+    Vector const half1 = rk4_step(f, position, static_cast<Scalar>(0.5) * step_size);
+    Vector const half2 = rk4_step(f, half1, static_cast<Scalar>(0.5) * step_size);
+
+    // TODO: I'm hard coding in an assumption that a Vector is only necessary because we're solving
+    // a second (or higher) order system, and so the only error we really care about is the final
+    // value. This should be more generic.
+    Scalar const error_estimate = std::abs((full_step - half2)[0]);
+    Scalar next_step_size = step_size;
+    if (error_estimate > tolerance) {
+        next_step_size *= smaller;
+    }
+    else if (error_estimate < static_cast<Scalar>(0.5) * tolerance) {
+        next_step_size *= bigger;
+    }
+
+    return {half2, next_step_size};
+}
+
 #endif

_code/pset_04/integrators/main.cpp

-#include <iostream>
 #include <Eigen/Dense>
+#include <iostream>
+#include <tuple>
 #include <vector>
 #include "integrators.h"
 
@@ -59,6 +60,31 @@ IntegrationPath<Scalar, Vector> integration_test(
     return result;
 }
 
+//--------------------------------------------------------------------------------------------------
+template <typename Integrator>
+IntegrationPath<Scalar, Vector> adaptive_integration_test(
+    Integrator const& integrator,
+    Vector const& x_initial,
+    Scalar t_final,
+    Scalar step_size
+) {
+    Vector x = x_initial;
+    Scalar t = 0;
+
+    IntegrationPath<Scalar, Vector> result;
+    uint32_t const n_steps = static_cast<uint32_t>(t_final / step_size) + 1;
+    result.reserve(n_steps);
+    result.push_back(t, x);
+
+    while (t < t_final) {
+        t += step_size;
+        std::tie(x, step_size) = integrator(x, step_size);
+        result.push_back(t, x);
+    }
+
+    return result;
+}
+
 //--------------------------------------------------------------------------------------------------
 int main() {
     HarmonicOscillator oscillator;
@@ -72,8 +98,14 @@ int main() {
         return rk4_step(oscillator, x, s);
     };
 
+    auto const adaptive_rk4_integrator = [&](Vector const& x, Scalar s) {
+        return adaptive_rk4_step(oscillator, x, s, 1e-3);
+    };
+
     auto euler_results = integration_test(euler_integrator, Vector(1, 0), 314.159265359, 0.01);
-    auto rk4_results = integration_test(rk4_integrator, Vector(1, 0), 314.159265359, 0.01);
     std::cout << euler_results.positions.back() << '\n';
+    auto rk4_results = integration_test(rk4_integrator, Vector(1, 0), 314.159265359, 0.01);
     std::cout << rk4_results.positions.back() << '\n';
+    auto adaptive_rk4_results = adaptive_integration_test(adaptive_rk4_integrator, Vector(1, 0), 314.159265359, 0.01);
+    std::cout << adaptive_rk4_results.positions.back() << '\n';
 }