1 /*
2 [auto_generated]
3 libs/numeric/odeint/test/bulirsch_stoer.cpp
4
5 [begin_description]
6 This file tests the Bulirsch-Stoer stepper.
7 [end_description]
8
9 Copyright 2011 Mario Mulansky
10 Copyright 2012 Karsten Ahnert
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 // disable checked iterator warning for msvc
19 #include <boost/config.hpp>
20 #ifdef BOOST_MSVC
21 #pragma warning(disable:4996)
22 #endif
23
24 #define BOOST_TEST_MODULE odeint_bulirsch_stoer
25
26 #include <utility>
27 #include <iostream>
28
29 #include <boost/array.hpp>
30
31 #include <boost/test/unit_test.hpp>
32
33 #include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
34 #include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
35
36 #include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
37
38 using namespace boost::unit_test;
39 using namespace boost::numeric::odeint;
40
41 typedef double value_type;
42 typedef boost::array< value_type , 3 > state_type;
43
44 const double sigma = 10.0;
45 const double R = 28.0;
46 const double b = 8.0 / 3.0;
47
48 struct lorenz
49 {
50 template< class State , class Deriv >
operator ()lorenz51 void operator()( const State &x , Deriv &dxdt , double t ) const
52 {
53 dxdt[0] = sigma * ( x[1] - x[0] );
54 dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
55 dxdt[2] = -b * x[2] + x[0] * x[1];
56 }
57 };
58
59 struct const_system
60 {
61 template< class State , class Deriv >
operator ()const_system62 void operator()( const State &x , Deriv &dxdt , double t ) const
63 {
64 dxdt[0] = 1.0;
65 dxdt[1] = 1.0;
66 dxdt[2] = 1.0;
67 }
68 };
69
70 struct sin_system
71 {
72 template< class State , class Deriv >
operator ()sin_system73 void operator()( const State &x , Deriv &dxdt , double t ) const
74 {
75 dxdt[0] = sin( x[0] );
76 dxdt[1] = cos( x[1] );
77 dxdt[2] = sin( x[2] ) + cos( x[2] );
78 }
79 };
80
81 BOOST_AUTO_TEST_SUITE( bulirsch_stoer_test )
82
BOOST_AUTO_TEST_CASE(test_bulirsch_stoer)83 BOOST_AUTO_TEST_CASE( test_bulirsch_stoer )
84 {
85 typedef bulirsch_stoer< state_type > stepper_type;
86 stepper_type stepper( 1E-9 , 1E-9 , 1.0 , 0.0 );
87
88 state_type x;
89 x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
90
91 double dt = 0.1;
92
93 //stepper.try_step( lorenz() , x , t , dt );
94
95 std::cout << "starting integration..." << std::endl;
96
97 size_t steps = integrate_adaptive( stepper , lorenz() , x , 0.0 , 10.0 , dt );
98
99 std::cout << "required steps: " << steps << std::endl;
100
101 bulirsch_stoer_dense_out< state_type > bs_do( 1E-9 , 1E-9 , 1.0 , 0.0 );
102 x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
103 double t = 0.0;
104 dt = 1E-1;
105 bs_do.initialize( x , t , dt );
106 bs_do.do_step( sin_system() );
107 std::cout << "one step successful, new time: " << bs_do.current_time() << " (" << t << ")" << std::endl;
108
109 x = bs_do.current_state();
110 std::cout << "x( " << bs_do.current_time() << " ) = [ " << x[0] << " , " << x[1] << " , " << x[2] << " ]" << std::endl;
111
112 bs_do.calc_state( bs_do.current_time()/3 , x );
113 std::cout << "x( " << bs_do.current_time()/3 << " ) = [ " << x[0] << " , " << x[1] << " , " << x[2] << " ]" << std::endl;
114
115 std::cout << std::endl << "=======================================================================" << std::endl << std::endl;
116
117 x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
118 t = 0.0; dt /= 3;
119 bs_do.initialize( x , t , dt );
120 bs_do.do_step( sin_system() );
121 x = bs_do.current_state();
122 std::cout << "x( " << bs_do.current_time() << " ) = [ " << x[0] << " , " << x[1] << " , " << x[2] << " ]" << std::endl;
123
124 t = dt;
125 bs_do.initialize( x , t , dt );
126 bs_do.do_step( sin_system() );
127 x = bs_do.current_state();
128
129 t = 2*dt;
130 bs_do.initialize( x , t , dt );
131 bs_do.do_step( sin_system() );
132 x = bs_do.current_state();
133
134 std::cout << "x( " << bs_do.current_time() << " ) = [ " << x[0] << " , " << x[1] << " , " << x[2] << " ]" << std::endl << std::endl << std::endl;
135 }
136
BOOST_AUTO_TEST_CASE(test_bulirsch_stoer_adjust_size)137 BOOST_AUTO_TEST_CASE( test_bulirsch_stoer_adjust_size )
138 {
139 typedef bulirsch_stoer< state_type > stepper_type;
140 stepper_type stepper( 1E-9 , 1E-9 , 1.0 , 0.0 );
141
142 state_type x; x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
143
144 stepper.adjust_size( x );
145
146
147 double dt = 0.1;
148
149 size_t steps = integrate_adaptive( stepper , lorenz() , x , 0.0 , 10.0 , dt );
150
151 std::cout << "required steps: " << steps << std::endl;
152 }
153
154
155 BOOST_AUTO_TEST_SUITE_END()
156