1 /*
2 * find_crossing.cpp
3 *
4 * Finds the energy threshold crossing for a damped oscillator.
5 * The algorithm uses a dense out stepper with find_if to first find an
6 * interval containing the threshold crossing and the utilizes the dense out
7 * functionality with a bisection to further refine the interval until some
8 * desired precision is reached.
9 *
10 * Copyright 2015 Mario Mulansky
11 *
12 * Distributed under the Boost Software License, Version 1.0.
13 * (See accompanying file LICENSE_1_0.txt or
14 * copy at http://www.boost.org/LICENSE_1_0.txt)
15 */
16
17
18
19 #include <iostream>
20 #include <utility>
21 #include <algorithm>
22 #include <array>
23
24 #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
25 #include <boost/numeric/odeint/stepper/generation.hpp>
26 #include <boost/numeric/odeint/iterator/adaptive_iterator.hpp>
27
28 namespace odeint = boost::numeric::odeint;
29
30 typedef std::array<double, 2> state_type;
31
32 const double gam = 1.0; // damping strength
33
damped_osc(const state_type & x,state_type & dxdt,const double)34 void damped_osc(const state_type &x, state_type &dxdt, const double /*t*/)
35 {
36 dxdt[0] = x[1];
37 dxdt[1] = -x[0] - gam * x[1];
38 }
39
40
41 struct energy_condition {
42
43 // defines the threshold crossing in terms of a boolean functor
44
45 double m_min_energy;
46
energy_conditionenergy_condition47 energy_condition(const double min_energy)
48 : m_min_energy(min_energy) { }
49
energyenergy_condition50 double energy(const state_type &x) {
51 return 0.5 * x[1] * x[1] + 0.5 * x[0] * x[0];
52 }
53
operator ()energy_condition54 bool operator()(const state_type &x) {
55 // becomes true if the energy becomes smaller than the threshold
56 return energy(x) <= m_min_energy;
57 }
58 };
59
60
61 template<class System, class Condition>
62 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)63 find_condition(state_type &x0, System sys, Condition cond,
64 const double t_start, const double t_end, const double dt,
65 const double precision = 1E-6) {
66
67 // integrates an ODE until some threshold is crossed
68 // returns time and state at the point of the threshold crossing
69 // if no threshold crossing is found, some time > t_end is returned
70
71 auto stepper = odeint::make_dense_output(1.0e-6, 1.0e-6,
72 odeint::runge_kutta_dopri5<state_type>());
73
74 auto ode_range = odeint::make_adaptive_range(std::ref(stepper), sys, x0,
75 t_start, t_end, dt);
76
77 // find the step where the condition changes
78 auto found_iter = std::find_if(ode_range.first, ode_range.second, cond);
79
80 if(found_iter == ode_range.second)
81 {
82 // no threshold crossing -> return time after t_end and ic
83 return std::make_pair(t_end + dt, x0);
84 }
85
86 // the dense out stepper now covers the interval where the condition changes
87 // improve the solution by bisection
88 double t0 = stepper.previous_time();
89 double t1 = stepper.current_time();
90 double t_m;
91 state_type x_m;
92 // use odeint's resizing functionality to allocate memory for x_m
93 odeint::adjust_size_by_resizeability(x_m, x0,
94 typename odeint::is_resizeable<state_type>::type());
95 while(std::abs(t1 - t0) > precision) {
96 t_m = 0.5 * (t0 + t1); // get the mid point time
97 stepper.calc_state(t_m, x_m); // obtain the corresponding state
98 if (cond(x_m))
99 t1 = t_m; // condition changer lies before midpoint
100 else
101 t0 = t_m; // condition changer lies after midpoint
102 }
103 // we found the interval of size eps, take it's midpoint as final guess
104 t_m = 0.5 * (t0 + t1);
105 stepper.calc_state(t_m, x_m);
106 return std::make_pair(t_m, x_m);
107 }
108
109
main(int argc,char ** argv)110 int main(int argc, char **argv)
111 {
112 state_type x0 = {{10.0, 0.0}};
113 const double t_start = 0.0;
114 const double t_end = 10.0;
115 const double dt = 0.1;
116 const double threshold = 0.1;
117
118 energy_condition cond(threshold);
119 state_type x_cond;
120 double t_cond;
121 std::tie(t_cond, x_cond) = find_condition(x0, damped_osc, cond,
122 t_start, t_end, dt, 1E-6);
123 if(t_cond > t_end)
124 {
125 // time after t_end -> no threshold crossing within [t_start, t_end]
126 std::cout << "No threshold crossing found." << std::endl;
127 } else
128 {
129 std::cout.precision(16);
130 std::cout << "Time of energy threshold crossing: " << t_cond << std::endl;
131 std::cout << "State: [" << x_cond[0] << " , " << x_cond[1] << "]" << std::endl;
132 std::cout << "Energy: " << cond.energy(x_cond) << std::endl;
133 }
134 }
135