1 /*
2 [auto_generated]
3 libs/numeric/odeint/examples/heun.cpp
4
5 [begin_description]
6 Examplary implementation of the method of Heun.
7 [end_description]
8
9 Copyright 2012 Karsten Ahnert
10 Copyright 2012 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 #include <iostream>
18
19
20 #include <boost/fusion/container/vector.hpp>
21 #include <boost/fusion/container/generation/make_vector.hpp>
22
23 #include <boost/array.hpp>
24
25 #include <boost/numeric/odeint.hpp>
26
27
28
29
30
31
32 namespace fusion = boost::fusion;
33
34 //[ heun_define_coefficients
35 template< class Value = double >
36 struct heun_a1 : boost::array< Value , 1 > {
heun_a1heun_a137 heun_a1( void )
38 {
39 (*this)[0] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
40 }
41 };
42
43 template< class Value = double >
44 struct heun_a2 : boost::array< Value , 2 >
45 {
heun_a2heun_a246 heun_a2( void )
47 {
48 (*this)[0] = static_cast< Value >( 0 );
49 (*this)[1] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
50 }
51 };
52
53
54 template< class Value = double >
55 struct heun_b : boost::array< Value , 3 >
56 {
heun_bheun_b57 heun_b( void )
58 {
59 (*this)[0] = static_cast<Value>( 1 ) / static_cast<Value>( 4 );
60 (*this)[1] = static_cast<Value>( 0 );
61 (*this)[2] = static_cast<Value>( 3 ) / static_cast<Value>( 4 );
62 }
63 };
64
65 template< class Value = double >
66 struct heun_c : boost::array< Value , 3 >
67 {
heun_cheun_c68 heun_c( void )
69 {
70 (*this)[0] = static_cast< Value >( 0 );
71 (*this)[1] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
72 (*this)[2] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
73 }
74 };
75 //]
76
77
78 //[ heun_stepper_definition
79 template<
80 class State ,
81 class Value = double ,
82 class Deriv = State ,
83 class Time = Value ,
84 class Algebra = boost::numeric::odeint::range_algebra ,
85 class Operations = boost::numeric::odeint::default_operations ,
86 class Resizer = boost::numeric::odeint::initially_resizer
87 >
88 class heun : public
89 boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
90 Algebra , Operations , Resizer >
91 {
92
93 public:
94
95 typedef boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
96 Algebra , Operations , Resizer > stepper_base_type;
97
98 typedef typename stepper_base_type::state_type state_type;
99 typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
100 typedef typename stepper_base_type::value_type value_type;
101 typedef typename stepper_base_type::deriv_type deriv_type;
102 typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
103 typedef typename stepper_base_type::time_type time_type;
104 typedef typename stepper_base_type::algebra_type algebra_type;
105 typedef typename stepper_base_type::operations_type operations_type;
106 typedef typename stepper_base_type::resizer_type resizer_type;
107 typedef typename stepper_base_type::stepper_type stepper_type;
108
heun(const algebra_type & algebra=algebra_type ())109 heun( const algebra_type &algebra = algebra_type() )
110 : stepper_base_type(
111 fusion::make_vector(
112 heun_a1<Value>() ,
113 heun_a2<Value>() ) ,
114 heun_b<Value>() , heun_c<Value>() , algebra )
115 { }
116 };
117 //]
118
119
120 const double sigma = 10.0;
121 const double R = 28.0;
122 const double b = 8.0 / 3.0;
123
124 struct lorenz
125 {
126 template< class State , class Deriv >
operator ()lorenz127 void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
128 {
129 typename boost::range_iterator< const State >::type x = boost::begin( x_ );
130 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
131
132 dxdt[0] = sigma * ( x[1] - x[0] );
133 dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
134 dxdt[2] = -b * x[2] + x[0] * x[1];
135 }
136 };
137
138 struct streaming_observer
139 {
140 std::ostream &m_out;
streaming_observerstreaming_observer141 streaming_observer( std::ostream &out ) : m_out( out ) { }
142 template< typename State , typename Value >
operator ()streaming_observer143 void operator()( const State &x , Value t ) const
144 {
145 m_out << t;
146 for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i];
147 m_out << "\n";
148 }
149 };
150
151
152
main(int argc,char ** argv)153 int main( int argc , char **argv )
154 {
155 using namespace std;
156 using namespace boost::numeric::odeint;
157
158
159 //[ heun_example
160 typedef boost::array< double , 3 > state_type;
161 heun< state_type > h;
162 state_type x = {{ 10.0 , 10.0 , 10.0 }};
163
164 integrate_const( h , lorenz() , x , 0.0 , 100.0 , 0.01 ,
165 streaming_observer( std::cout ) );
166
167 //]
168
169 return 0;
170 }
171