1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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;
}
}
|
