1 /*
2 * Simulation of an ensemble of Roessler attractors
3 *
4 * Copyright 2014 Mario Mulansky
5 *
6 * Distributed under the Boost Software License, Version 1.0.
7 * (See accompanying file LICENSE_1_0.txt or
8 * copy at http://www.boost.org/LICENSE_1_0.txt)
9 *
10 */
11
12
13 #include <iostream>
14 #include <vector>
15 #include <random>
16
17 #include <boost/timer.hpp>
18 #include <boost/array.hpp>
19
20 #include <boost/numeric/odeint.hpp>
21
22 namespace odeint = boost::numeric::odeint;
23
24 typedef boost::timer timer_type;
25
26 typedef double fp_type;
27 //typedef float fp_type;
28
29 typedef boost::array<fp_type, 3> state_type;
30 typedef std::vector<state_type> state_vec;
31
32 //---------------------------------------------------------------------------
33 struct roessler_system {
34 const fp_type m_a, m_b, m_c;
35
roessler_systemroessler_system36 roessler_system(const fp_type a, const fp_type b, const fp_type c)
37 : m_a(a), m_b(b), m_c(c)
38 {}
39
operator ()roessler_system40 void operator()(const state_type &x, state_type &dxdt, const fp_type t) const
41 {
42 dxdt[0] = -x[1] - x[2];
43 dxdt[1] = x[0] + m_a * x[1];
44 dxdt[2] = m_b + x[2] * (x[0] - m_c);
45 }
46 };
47
48 //---------------------------------------------------------------------------
main(int argc,char * argv[])49 int main(int argc, char *argv[]) {
50 if(argc<3)
51 {
52 std::cerr << "Expected size and steps as parameter" << std::endl;
53 exit(1);
54 }
55 const size_t n = atoi(argv[1]);
56 const size_t steps = atoi(argv[2]);
57 //const size_t steps = 50;
58
59 const fp_type dt = 0.01;
60
61 const fp_type a = 0.2;
62 const fp_type b = 1.0;
63 const fp_type c = 9.0;
64
65 // random initial conditions on the device
66 std::vector<fp_type> x(n), y(n), z(n);
67 std::default_random_engine generator;
68 std::uniform_real_distribution<fp_type> distribution_xy(-8.0, 8.0);
69 std::uniform_real_distribution<fp_type> distribution_z(0.0, 20.0);
70 auto rand_xy = std::bind(distribution_xy, std::ref(generator));
71 auto rand_z = std::bind(distribution_z, std::ref(generator));
72 std::generate(x.begin(), x.end(), rand_xy);
73 std::generate(y.begin(), y.end(), rand_xy);
74 std::generate(z.begin(), z.end(), rand_z);
75
76 state_vec state(n);
77 for(size_t i=0; i<n; ++i)
78 {
79 state[i][0] = x[i];
80 state[i][1] = y[i];
81 state[i][2] = z[i];
82 }
83
84 std::cout.precision(16);
85
86 std::cout << "# n: " << n << std::endl;
87
88 std::cout << x[0] << std::endl;
89
90
91 // Stepper type - use never_resizer for slight performance improvement
92 odeint::runge_kutta4_classic<state_type, fp_type, state_type, fp_type,
93 odeint::array_algebra,
94 odeint::default_operations,
95 odeint::never_resizer> stepper;
96
97 roessler_system sys(a, b, c);
98
99 timer_type timer;
100
101 fp_type t = 0.0;
102
103 for (int step = 0; step < steps; step++)
104 {
105 for(size_t i=0; i<n; ++i)
106 {
107 stepper.do_step(sys, state[i], t, dt);
108 }
109 t += dt;
110 }
111
112 std::cout << "Integration finished, runtime for " << steps << " steps: ";
113 std::cout << timer.elapsed() << " s" << std::endl;
114
115 // compute some accumulation to make sure all results have been computed
116 fp_type s = 0.0;
117 for(size_t i = 0; i < n; ++i)
118 {
119 s += state[i][0];
120 }
121
122 std::cout << state[0][0] << std::endl;
123 std::cout << s/n << std::endl;
124
125 }
126