1 /*
2 * odeint_rk4_array
3 *
4 * Copyright 2011 Mario Mulansky
5 * Copyright 2012 Karsten Ahnert
6 *
7 * Distributed under the Boost Software License, Version 1.0.
8 * (See accompanying file LICENSE_1_0.txt or
9 * copy at http://www.boost.org/LICENSE_1_0.txt)
10 */
11
12 #include <iostream>
13
14 #include <boost/timer.hpp>
15 #include <boost/array.hpp>
16
17 #include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
18 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
19 #include <boost/numeric/odeint/algebra/array_algebra.hpp>
20
21 #include "lorenz.hpp"
22
23 typedef boost::timer timer_type;
24
25 typedef boost::array< double , 3 > state_type;
26
27 using namespace boost::numeric::odeint;
28
29 //typedef boost::numeric::odeint::runge_kutta4_classic< state_type > rk4_odeint_type;
30
31 // use the never resizer explicitely for optimal performance with gcc,
32 // for the intel compiler this doesnt matter and the above definition
33 // gives the same performance
34 typedef runge_kutta4_classic< state_type , double , state_type , double ,
35 array_algebra, default_operations, never_resizer > rk4_odeint_type;
36
37
38 const int loops = 21;
39 const int num_of_steps = 20000000;
40 const double dt = 1E-10;
41
42
main()43 int main()
44 {
45 double min_time = 1E6; // something big
46 rk4_odeint_type stepper;
47 std::clog.precision(16);
48 std::cout.precision(16);
49 for( int n=0; n<loops; n++ )
50 {
51 state_type x = {{ 8.5, 3.1, 1.2 }};
52 double t = 0.0;
53 timer_type timer;
54 for( size_t i = 0 ; i < num_of_steps ; ++i )
55 {
56 stepper.do_step( lorenz(), x, t, dt );
57 t += dt;
58 }
59 min_time = std::min( timer.elapsed() , min_time );
60 std::clog << timer.elapsed() << '\t' << x[0] << std::endl;
61 }
62 std::cout << "Minimal Runtime: " << min_time << std::endl;
63 }
64