1 files changed, 134 insertions, 0 deletions
diff --git a/src/boost/libs/numeric/odeint/examples/find_crossing.cpp b/src/boost/libs/numeric/odeint/examples/find_crossing.cpp
new file mode 100644
index 00000000..a17ef4f9
--- /dev/null
+++ b/src/boost/libs/numeric/odeint/examples/find_crossing.cpp
@@ -0,0 +1,134 @@
+/*
+ * find_crossing.cpp
+ *
+ * Finds the energy threshold crossing for a damped oscillator.
+ * The algorithm uses a dense out stepper with find_if to first find an
+ * interval containing the threshold crossing and the utilizes the dense out
+ * functionality with a bisection to further refine the interval until some
+ * desired precision is reached.
+ *
+ * Copyright 2015 Mario Mulansky
+ *
+ * Distributed under the Boost Software License, Version 1.0.
+ * (See accompanying file LICENSE_1_0.txt or
+ * copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+
+
+#include <iostream>
+#include <utility>
+#include <algorithm>
+#include <array>
+
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/generation.hpp>
+#include <boost/numeric/odeint/iterator/adaptive_iterator.hpp>
+
+namespace odeint = boost::numeric::odeint;
+
+typedef std::array<double, 2> state_type;
+
+const double gam = 1.0;  // damping strength
+
+void damped_osc(const state_type &x, state_type &dxdt, const double /*t*/)
+{
+    dxdt[0] = x[1];
+    dxdt[1] = -x[0] - gam * x[1];
+}
+
+
+struct energy_condition {
+
+    // defines the threshold crossing in terms of a boolean functor
+
+    double m_min_energy;
+
+    energy_condition(const double min_energy)
+        : m_min_energy(min_energy) { }
+
+    double energy(const state_type &x) {
+        return 0.5 * x[1] * x[1] + 0.5 * x[0] * x[0];
+    }
+
+    bool operator()(const state_type &x) {
+        // becomes true if the energy becomes smaller than the threshold
+        return energy(x) <= m_min_energy;
+    }
+};
+
+
+template<class System, class Condition>
+std::pair<double, state_type>
+find_condition(state_type &x0, System sys, Condition cond,
+               const double t_start, const double t_end, const double dt,
+               const double precision = 1E-6) {
+
+    // integrates an ODE until some threshold is crossed
+    // returns time and state at the point of the threshold crossing
+    // if no threshold crossing is found, some time > t_end is returned
+
+    auto stepper = odeint::make_dense_output(1.0e-6, 1.0e-6,
+                                             odeint::runge_kutta_dopri5<state_type>());
+
+    auto ode_range = odeint::make_adaptive_range(std::ref(stepper), sys, x0,
+                                                 t_start, t_end, dt);
+
+    // find the step where the condition changes
+    auto found_iter = std::find_if(ode_range.first, ode_range.second, cond);
+
+    if(found_iter == ode_range.second)
+    {
+        // no threshold crossing -> return time after t_end and ic
+        return std::make_pair(t_end + dt, x0);
+    }
+
+    // the dense out stepper now covers the interval where the condition changes
+    // improve the solution by bisection
+    double t0 = stepper.previous_time();
+    double t1 = stepper.current_time();
+    double t_m;
+    state_type x_m;
+    // use odeint's resizing functionality to allocate memory for x_m
+    odeint::adjust_size_by_resizeability(x_m, x0,
+                                         typename odeint::is_resizeable<state_type>::type());
+    while(std::abs(t1 - t0) > precision) {
+        t_m = 0.5 * (t0 + t1);  // get the mid point time
+        stepper.calc_state(t_m, x_m); // obtain the corresponding state
+        if (cond(x_m))
+            t1 = t_m;  // condition changer lies before midpoint
+        else
+            t0 = t_m;  // condition changer lies after midpoint
+    }
+    // we found the interval of size eps, take it's midpoint as final guess
+    t_m = 0.5 * (t0 + t1);
+    stepper.calc_state(t_m, x_m);
+    return std::make_pair(t_m, x_m);
+}
+
+
+int main(int argc, char **argv)
+{
+    state_type x0 = {{10.0, 0.0}};
+    const double t_start = 0.0;
+    const double t_end = 10.0;
+    const double dt = 0.1;
+    const double threshold = 0.1;
+
+    energy_condition cond(threshold);
+    state_type x_cond;
+    double t_cond;
+    std::tie(t_cond, x_cond) = find_condition(x0, damped_osc, cond,
+                                              t_start, t_end, dt, 1E-6);
+    if(t_cond > t_end)
+    {
+        // time after t_end -> no threshold crossing within [t_start, t_end]
+        std::cout << "No threshold crossing found." << std::endl;
+    } else
+    {
+        std::cout.precision(16);
+        std::cout << "Time of energy threshold crossing: " << t_cond << std::endl;
+        std::cout << "State: [" << x_cond[0] << " , " << x_cond[1] << "]" << std::endl;
+        std::cout << "Energy: " << cond.energy(x_cond) << std::endl;
+    }
+}