1 /*
2 * stepper_details.cpp
3 *
4 * This example demonstrates some details about the steppers in odeint.
5 *
6 *
7 * Copyright 2011-2012 Karsten Ahnert
8 * Copyright 2012 Mario Mulansky
9 * Copyright 2013 Pascal Germroth
10 *
11 * Distributed under the Boost Software License, Version 1.0.
12 * (See accompanying file LICENSE_1_0.txt or
13 * copy at http://www.boost.org/LICENSE_1_0.txt)
14 */
15
16 #include <iostream>
17 #include <boost/array.hpp>
18 #include <boost/bind.hpp>
19 #include <boost/numeric/odeint.hpp>
20
21 using namespace std;
22 using namespace boost::numeric::odeint;
23
24 const size_t N = 3;
25
26 typedef boost::array< double , N > state_type;
27
28 //[ system_function_structure
sys(const state_type &,state_type &,const double)29 void sys( const state_type & /*x*/ , state_type & /*dxdt*/ , const double /*t*/ )
30 {
31 // ...
32 }
33 //]
34
sys1(const state_type &,state_type &,const double)35 void sys1( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
36 {
37 }
38
sys2(const state_type &,state_type &,const double)39 void sys2( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
40 {
41 }
42
43
44 //[ symplectic_stepper_detail_system_function
45 typedef boost::array< double , 1 > vector_type;
46
47
48 struct harm_osc_f1
49 {
operator ()harm_osc_f150 void operator()( const vector_type &p , vector_type &dqdt )
51 {
52 dqdt[0] = p[0];
53 }
54 };
55
56 struct harm_osc_f2
57 {
operator ()harm_osc_f258 void operator()( const vector_type &q , vector_type &dpdt )
59 {
60 dpdt[0] = -q[0];
61 }
62 };
63 //]
64
65 //[ symplectic_stepper_detail_system_class
66 struct harm_osc
67 {
f1harm_osc68 void f1( const vector_type &p , vector_type &dqdt ) const
69 {
70 dqdt[0] = p[0];
71 }
72
f2harm_osc73 void f2( const vector_type &q , vector_type &dpdt ) const
74 {
75 dpdt[0] = -q[0];
76 }
77 };
78 //]
79
main(int argc,char ** argv)80 int main( int argc , char **argv )
81 {
82 using namespace std;
83
84 // Explicit stepper example
85 {
86 double t( 0.0 ) , dt( 0.1 );
87 state_type in , out , dxdtin , inout;
88 //[ explicit_stepper_detail_example
89 runge_kutta4< state_type > rk;
90 rk.do_step( sys1 , inout , t , dt ); // In-place transformation of inout
91 rk.do_step( sys2 , inout , t , dt ); // call with different system: Ok
92 rk.do_step( sys1 , in , t , out , dt ); // Out-of-place transformation
93 rk.do_step( sys1 , inout , dxdtin , t , dt ); // In-place tranformation of inout
94 rk.do_step( sys1 , in , dxdtin , t , out , dt ); // Out-of-place transformation
95 //]
96 }
97
98
99
100 // FSAL stepper example
101 {
102 double t( 0.0 ) , dt( 0.1 );
103 state_type in , in2 , in3 , out , dxdtin , dxdtout , inout , dxdtinout;
104 //[ fsal_stepper_detail_example
105 runge_kutta_dopri5< state_type > rk;
106 rk.do_step( sys1 , in , t , out , dt );
107 rk.do_step( sys2 , in , t , out , dt ); // DONT do this, sys1 is assumed
108
109 rk.do_step( sys2 , in2 , t , out , dt );
110 rk.do_step( sys2 , in3 , t , out , dt ); // DONT do this, in2 is assumed
111
112 rk.do_step( sys1 , inout , dxdtinout , t , dt );
113 rk.do_step( sys2 , inout , dxdtinout , t , dt ); // Ok, internal derivative is not used, dxdtinout is updated
114
115 rk.do_step( sys1 , in , dxdtin , t , out , dxdtout , dt );
116 rk.do_step( sys2 , in , dxdtin , t , out , dxdtout , dt ); // Ok, internal derivative is not used
117 //]
118 }
119
120
121 // Symplectic harmonic oscillator example
122 {
123 double t( 0.0 ) , dt( 0.1 );
124 //[ symplectic_stepper_detail_example
125 pair< vector_type , vector_type > x;
126 x.first[0] = 1.0; x.second[0] = 0.0;
127 symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
128 rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , x , t , dt );
129 //]
130
131 //[ symplectic_stepper_detail_system_class_example
132 harm_osc h;
133 rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) ,
134 x , t , dt );
135 //]
136 }
137
138 // Simplified harmonic oscillator example
139 {
140 double t = 0.0, dt = 0.1;
141 //[ simplified_symplectic_stepper_example
142 pair< vector_type , vector_type > x;
143 x.first[0] = 1.0; x.second[0] = 0.0;
144 symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
145 rkn.do_step( harm_osc_f1() , x , t , dt );
146 //]
147
148 vector_type q = {{ 1.0 }} , p = {{ 0.0 }};
149 //[ symplectic_stepper_detail_ref_usage
150 rkn.do_step( harm_osc_f1() , make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt );
151 rkn.do_step( harm_osc_f1() , q , p , t , dt );
152 rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , q , p , t , dt );
153 //]
154 }
155
156 // adams_bashforth_moulton stepper example
157 {
158 double t = 0.0 , dt = 0.1;
159 state_type inout;
160 //[ multistep_detail_example
161 adams_bashforth_moulton< 5 , state_type > abm;
162 abm.initialize( sys , inout , t , dt );
163 abm.do_step( sys , inout , t , dt );
164 //]
165
166 //[ multistep_detail_own_stepper_initialization
167 abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt );
168 //]
169 }
170
171
172
173 // dense output stepper examples
174 {
175 double t = 0.0 , dt = 0.1;
176 state_type in;
177 //[ dense_output_detail_example
178 dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense;
179 dense.initialize( in , t , dt );
180 pair< double , double > times = dense.do_step( sys );
181 (void)times;
182 //]
183
184 state_type inout;
185 double t_start = 0.0 , t_end = 1.0;
186 //[ dense_output_detail_generation1
187 typedef boost::numeric::odeint::result_of::make_dense_output<
188 runge_kutta_dopri5< state_type > >::type dense_stepper_type;
189 dense_stepper_type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
190 (void)dense2;
191 //]
192
193 //[ dense_output_detail_generation2
194 integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt );
195 //]
196 }
197
198
199 return 0;
200 }
201